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Foreword 


This  report  summarizes  work  done  on  AFGL  Contract  No.  F19628-81-1C-0046. 
While  specific  goals  of  the  research  are  discussed  in  Chapter  I,  it  is 
worthwhile  to  briefly  discuss  this  work  from  a  philosophical  and  strategical 
point  of  view. 

Firstly,  one  must  recognize  that  when  considering  an  entire  hemisphere  on  a 
routine  basis,  a  wide  variety  of  boundary  layer  situations  are  encountered.  In 
contrast,  boundary  layer  models  have  been  developed  for  special  situations  such 
as  the  convective  well-mixed  layer  or  the  fully  turbulent  stable  boundary 
layer.  In  the  real  world,  these  models  represent  asymptotic  limits  or  special 
cases . 

Secondly,  a  given  grid  box  contains  a  variety  of  surface  types  and 
associated  boundary  layer  responses.  Advection  from  one  part  of  the  grid  box  to 
the  other  is  important  in  the  corresponding  real  world.  In  the  model,  this 
advection  can  appear  only  as  subgrid  scale  flux.  Since  the  flux 
parameterization  is  necessarily  based  on  turbulence  theory,  such  advection  is  in 
effect  omitted  in  present  formulations. 

As  a  result  of  these  complications,  and  others,  only  a  simple  boundary 
layer  model  is  appropriate;  additional  sophistication  is  wasted  no  matter  how 
appealing  from  a  physical  point  of  view.  In  addition  to  simplicity,  the  model 
must  be  robust.  It  cannot  breakdown  or  lead  to  instability  in  special 
situations  even  if  rare.  Thus  we  are  willing  to  tolerate  substantial  errors  or 
physical  shortcomings  in  order  to  achieve  versatility.  This  is  not  the  usual 
initial  goal  of  modelling,  but  in  effect  becomes  the  goal  as  the  model  is  tested 
or  runs  operationally. 


In  the  development  of  our  model,  we  have  tried  to  keep  this  perspective  at 
the  forefront.  However,  our  model  contains  many  unique  features  which  we 
consider  to  be  important  improvements,  but  have  not  been  tested  under  a  variety 
of  conditions  nor  tested  on  a  routine  basis.  With  complete  certainty, 
troubleshooting  and  model  modification  will  be  necessary  after  use  with  the 
entire  atmospheric  model  under  a  variety  of  conditions. 
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I.  Introduction 


This  report  presents  the  physical  motivation  for  the  components  of  the  one 
dimensional  model  and  discusses  the  advantages  and  disadvantages  of  each 
component . 

The  main  components  of  the  model  are  the  evaporation  routine  (Chapter  II), 
the  soil  hydrology  model  (Chapter  III),  the  canopy-transpiration  model  (Chapter 
IV)  and  the  atmospheric  boundary  layer  model  (Chapter  V).  The  need  for 
simplicity  is  imposed  on  all  facets  of  the  model  development. 

The  organization  of  this  report  is  conveniently  illustrated  by  considering 
the  expression  for  total  evapotranspiration  of  water  from  the  soil-vegetation 
complex  to  the  atmosphere 

E  *  E. .  +E+E  *  P  •  E 
dir  T  c  p 

where  ^jir  is  the  direct  evaporation  from  the  soil,  is  the  transpiration, 

Ec  is  the  direct  evaporation  of  precipitation  intercepted  by  the  canopy  and 
Ep  is  the  potential  evaporation.  The  three  contributions  can  be  linearly 
related  to  the  potential  evaporation  and  functions  of  the  soil  water  distribu¬ 
tion,  plant  density,  stomatal  resistance  and  canopy  water  content.  F  is  the  sum 
of  such  functions. 

Chapter  II  discusses  the  representation  of  potential  evaporation  and 
devotes  considerable  attention  to  the  inclusion  of  the  influence  of  atmospheric 
stability  in  the  Penman  relationship  for  potential  evaporation.  The  stability 
influence  has  not  been  previously  included  at  the  level  of  simplicity  sought  in 
this  development.  The  use  of  24-hour  averaged  variables  is  examined  in  some 
detail  since  transpiration  parameters  are  typically  defined  in  terms  of  24-hour 
averaged  variables.  Over  water  surfaces,  we  are  recommending  the  usual  bulk 
aerodynamic  relationship  (Chapter  V). 


Oh;ipt-<»r  TIT  develops  the  two-layer  soil  hydrology  model  with  special 
emphasis  on  the  representation  of  transport  of  water  near  the  air-soil  interface 
and  on  the  behavior  of  truncation  errors  which  become  important  in  a  two-layer 
representation.  The  behavior  of  truncation  errors  is  studied  by  comparing 
results  with  a  high  resolution  model. 

The  soil  heat  flux  is  not  modelled  here.  Instead  a  surface  energy  balance 
is  employed  (Chapter  V) .  The  heat  flux  from  the  soil  to  the  atmosphere  is 
small,  particularly  when  accumulated  over  a  day  or  more.  Heat  transport  from 
the  soil  is  most  unimportant  with  organic  debris  or  litter  such  as  dead  grass, 
leaves  or  needles.  The  conductivity  of  such  a  layer  is  small,  especially  when 
dry.  From  a  global  point  of  view,  land  surfaces  are  typically  covered  with  such 
debris.  The  commonly-studied  bare  soil  case  occupies  only  a  small  percentage  of 
the  land  surface.  However,  even  with  organic  debris,  the  evapotranspiration  to 
the  atmosphere  may  still  be  large.  As  a  result,  water  flux  from  the  soil  to  the 
atmosphere  is  often  important  even  when  the  corresponding  heat  flux  is  not 
important. 

Chapter  IV  develops  the  transpiration  and  canopy  water  budget  relation¬ 
ships.  Section  IV. 1  motivates  the  need  for  a  canopy  water  budget  and  surveys 
various  approaches.  Section  IV. 2  discusses  the  representation  of  transpira¬ 
tion.  Section  IV. 3  presents  the  model  of  transpiration  and  canopy  water  budget 
used  here.  Chapter  V  develops  the  model  for  the  atmospheric  boundary  layer.  A 
new  depth  formulation  is  based  on  a  bulk  Richardson  number  which  is  generalized 
to  include  the  free  convection  case.  This  generalization  is  developed  concur¬ 
rently  with  the  countergradient  heat  flux  modification  to  the  eddy  exchange 
coefficient. 

The  values  of  transpiration  coefficients  and  soil  properties  appearing  in 
the  model  have  been  chosen  in  an  attempt  to  represent  typical  or  average  condi- 


tions  in  a  global  sense.  This  selection  is  tenuous  and  there  undoubtedly  exists 
a  prejudice  toward  mid-latitude  situations  where  most  of  the  observational 
evidence  has  been  collected.  Specification  of  the  global  distribution  of  such 
values  would  be  useful  future  work. 

He  emphasize  that  much  of  the  material  presented  in  this  report  motivates, 
develops,  and  tests  various  model  components.  For  purposes  of  quickly 
identifying  the  main  equations,  one  need  only  to  read  Sections  II. 2-3,  III. 2, 
III. 4,  IV. 3  and  V.  The  components  of  the  model  are  discussed  independently  so 
that  it  is  not  necessary  to  study  the  chapters  in  order.  The  software  and  users 
guide  for  operation  of  the  model  have  been  provided  to  the  Atmospheric 
Prediction  Branch  of  the  Air  Force  Geophysics  Laboratory. 


II.  The  Influence  of  Atmospheric  Stability  on  Potential  Evaporation 

Abstract 

The  Penman  relationship  for  potential  evaporation  is  modified  to  simply 
include  the  influence  of  atmospheric  stability  on  turbulent  transport  of  water 
vapor.  Explicit  expressions  for  the  stability-dependent,  surface  exchange 
coefficient  developed  by  Louis  are  used.  The  diurnal  variation  of  potential 
evaporation  is  computed  for  the  stability-dependent  and  original  Penman 
relationships  using  Wangara  data. 

The  influence  of  afternoon  instability  increases  the  aerodynamic  term  of 
the  modified  Penman  relationship  by  50%  or  more  on  days  with  moderate 
instability.  However,  the  unmodified  Penman  relationship  predicts  values  of 
daily  potential  evaporation  close  to  that  of  the  stability-dependent 
relationship.  This  agreement  is  partly  due  to  compensating  overestimation 
during  nighttime  hours.  Errors  due  to  use  of  daily-averaged  variables  are 
examined  in  detail  by  evaluating  the  nonlinear  interactions  between  the  diurnal 
variation  of  the  variables  in  the  Penman  relationship. 

A  simpler  method  for  estimating  the  exchange  coefficient  is  constructed 
from  an  empirical  relationship  between  the  radiation-Richardson  number  and  the 
Obukhov  length.  This  method  is  less  accurate  but  allows  estimation  of  the 
stability-dependent  exchange  coefficient  using  only  parameters  already  required 
for  evaluation  of  the  Penman  relationship.  Finally,  the  diurnal  variation  of 
the  atmospheric  resistance  coefficient  appearing  in  the  Penman-Monte ith 
relationship  is  presented. 


1.  Introduction 


The  surface  moisture  flux  is  often  parameterized  in  terms  of  potential 
evaporation  and  associated  coefficient  representing  soil  moisture  deficit, 
resistance  of  the  vegetation,  and  radiative  properties  of  the  surface.  The 
potential  evaporation  is  defined  as  that  evaporation  occurring  over  a  free  water 
surface.  In  theory,  the  potential  evaporation  is  independent  of  the  state  of 
the  water  surface  and  depends  only  on  atmospheric  conditions.  In  practice  the 
value  of  the  potential  evaporation  depends  on  the  methodology  of  measurement. 

In  most  atmospheric  models,  the  potential  evaporation  is  parameterized  in 
terms  of  bulk  aerodynamic  relationship  while  in  other  disciplines  the  potential 
evaporation  is  more  often  equated  to  a  Penman  (1948)  or  combination  formula¬ 
tion.  The  Penman  formulation  can  be  derived  by  combining  the  aerodynamic  rela¬ 
tionship  with  the  surface  energy  balance.  The  Penman  approach  has  the  following 
advantages: 

a)  Surface  temperature  or  surface  saturation  vapor  pressure  is  eliminat¬ 
ed.  In  practice  surface  temperature  is  difficult  to  define  over  land 
where  the  difference  between  vegetation  and  soil  temperatures  can 
exceed  5*C.  Associated  errors  in  the  bulk  aerodynamic  relationship 
have  been  shown  to  be  large  (Yu,  1977) . 

b)  The  Penman  relationship  includes  an  explicit  dependence  on  net  radia¬ 
tion  which,  when  calibrated  to  actual  evapotranspiration,  may  indirect¬ 
ly  include  biological  dependencies  on  solar  radiation  such  as  photosyn¬ 
thesis.  In  fact,  the  frequently  used  Priestley-Taylor  model  (1972) 
expresses  evapotranspiration  exclusively  in  terms  of  net  radiation. 

c)  The  Penman  relationship  has  been  compared  to  actual  evaporation  or 
evapotranspiration  in  a  large  number  of  studies,  although  many  studies 
are  hard  to  compare  due  to  use  of  different  versions  of  the  Penman 


equation,  use  of  different  observational  levels,  and  influences  of 
horizontal  inhomogeneity. 

In  comparison  with  more  empirical  approaches,  the  Penman  relationship  is 
usually  found  to  perform  as  well  or  better  (e.g.,  Seguin,  1975).  The  Penman 
relationship  has  not  been  tested  explicitly  for  downward  moisture  flux  associat¬ 
ed  with  negative  net  radiation  and  dew  formation. 

Evaluation  of  any  model  of  potential  evaporation  from  atmospheric  vari¬ 
ables,  which  are  in  equilibrium  with  a  surface  evaporating  at  less  than  poten¬ 
tial,  can  be  considered  to  be  inconsistent  (Bouchet,  1963;  Horton,  1975  and 
others).  Here  we  require  the  potential  evaporation  to  be  a  function  of  only 
atmospheric  variables  and  independent  of  reduction  of  actual  moisture  flux  due 
to  soil  moisture  deficit  and  plant  resistance.  However,  the  results  of  the 
present  development  can  be  transformed  to  modified  expressions  of  potential 
evaporation. 

The  most  serious  disadvantage  of  the  usual  Penman  relationship,  and  many 
other  models  of  potential  evaporation,  is  failure  to  explicitly  include  the 
influence  of  atmospheric  stability  cm  atmospheric  transport  of  water  vapor. 

Such  an  influence  can  significantly  contribute  to  diurnal  variation  of  the 
potential  evaporation.  The  influence  of  stability  can  be  reduced  by  using 
atmospheric  variables  measured  closer  to  the  ground.  However,  observations 
close  to  or  within  the  canopy  are  difficult  to  interpret  and  the  usual 
similarity  theory  no  longer  applies.  As  a  result,  studies  of  turbulent  fluxes 
over  land  almost  always  include  the  influence  of  atmospheric  stability. 

However,  only  a  few  of  the  many  applications  of  the  Penman  or  combination 
formulations  have  included  the  influence  of  atmospheric  stability.  Such 


formulations  are  usually  based  upon  similarity  modification  of  the  log  law  and 
thus  also  include  a  dependence  on  surface  roughness  in  contrast  to  the  original 
Penman  relationship.  Businger  (1956)  modified  the  Penman  relationship  to 
include  a  stability  correction  which  was  expressed  in  monogram  form.  Fuchs  et 
al.  (1969)  included  the  influence  of  stability  in  a  version  of  the  combination 
equation  which  did  not  eliminate  surface  temperature.  Federer  (1970)  included  a 
stability  adjustment  in  the  Penman  relationship  which  required  knowledge  of  the 
Obukhov  length  or  an  additional  unspecified  relationship  between  stability 
parameters . 

The  Penman-Monteith  (Monteith,  1965)  relationship  has  been  modified  to 
include  certain  aspects  of  the  stability  influence  using  Obukhov  similarity 
theory  (Stewart  and  Thom,  1973;  Thom  and  Oliver,  1977;  Verma  et  al.,  1976; 
DeHeer-Amiasah  et  al . ,  1981;  Berkowicz  and  Prahm,  1982).  Such  inclusion  of  the 
stability  influenoe  generally  requires  iterative  procedures.  Strieker  and 
Brutsaert  (1978)  apply  an  iterative  technique  to  estimate  the  stability  para¬ 
meter  and  its  influence  on  the  actual  surface  moisture  flux  as  computed  from  the 
surface  energy  budget,  while  Brutsaert  (1982)  suggests  an  iterative  procedure 
based  on  the  Penman  relationship.  They  conclude  that  the  influence  of  atmos¬ 
pheric  stability  cannot  be  neglected  although  they  found  little  difference 
between  the  various  stability  formulations  examined. 

The  first  goal  of  this  study  is  to  present  a  method  for  estimating  poten¬ 
tial  evaporation  from  the  Penman  relationship  which  includes  the  influence  of 
atmospheric  stability  yet  is  simpler  than  the  procedures  above.  In  particular, 
we  wish  to  avoid  the  need  for  iteration  in  order  to  construct  a  method  suitable 
for  use  in  those  atmospheric  models  or  routine  applications  where  computer  time 


is  restricted.  This  will  be  done  by  using  the  dependence  of  surface  exchange 
coefficients  on  the  bulk  Richardson  number  presented  in  Louis  (1979)  and  Louis 
et  al .  (1982)  for  both  the  stable  and  unstable  cases  (Section  2). 

A  second  goal  of  this  study  is  to  estimate  the  importance  of  the  influence 
of  atmospheric  stability  since  most  applications  of  the  Penman  relationship 
still  neglect  such  an  influence.  Toward  this  goal,  we  will  systematically 
evaluate  the  Penman  relationship  with  and  without  the  stability  influence  using 
data  from  the  Wangara  experiment  (Clarke  et  al.r  1971).  The  inclusion  of  the 
stability  influence  in  this  study  is  expected  to  substantially  improve  the 
Penman  relationships  since  the  Louis  (1979)  formulation  approximates  similarity 
theory  which  has  been  calibrated  and  tested  against  classical  data  sets.  How¬ 
ever,  the  modelled  stability  influence  will  incur  some  error  so  that  the  evalua¬ 
tion  of  the  original  Penman  relationship  will  be  only  approximate. 

Also  of  interest  is  the  influence  of  the  diurnal  variation  of  stability  on 
the  24-hour  evaporation  since  evaporation  values  are  often  reported  for  24-hour 
periods.  While  afternoon  instability  can  significantly  enhance  the  surface 
moisture  flux,  nocturnal  stability  can  significantly  reduce  such  fluxes  leading 
to  some  cancellation  between  stability  influences.  The  third  goal  of  this  study 
is  to  assess  the  magnitude  of  errors  resulting  from  use  of  24-hour  averaged 
variables  since  the  Penman  relationship  is  frequently  evaluated  with  such 
averaged  data. 

Note  that  this  study  is  concerned  only  with  potential  evaporation  dictated 
by  atmospheric  variables;  no  attempts  are  made  to  estimate  the  actual  evapora¬ 
tion  or  relate  it  to  the  potential  evaporation.  The  daily  potential  evaporation 
during  the  Wangara  observational  period  averaged  a  little  more  than  2mm  per  day 
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ceaching  near  4mn  on  some  days.  These  modest  values  would  be  more  than  enough 
to  evaporate  the  2  1/2  cm  of  precipitation  which  fell  during  the  43-day  period. 
Therefore,  the  actual  evaporation  rate  was  probably  well  below  the  potential 
rate  for  much  of  the  observational  period. 

The  study  of  the  behavior  of  potential  evaporation  under  conditions  where 
the  actual  evaporation  is  less  than  potential  is  of  considerable  interest. 
Plant-soil  models,  which  are  forced  by  expressions  for  potential  evaporation, 
are  typically  applied  to  situations  of  moisture  stress. 

It  should  be  noted  that  models  or  expressions  which  relate  actual  evapora¬ 
tion  to  potential  evaporation  depend  largely  on  observational  calibration. 
Consequently,  improved  physical  basis  sought  in  this  study  will  not  necessarily 
lead  to  improved  prediction  of  evapotranspiration  in  field  situations.  The 
relationship  between  actual  and  potential  evaporation  are  addressed  in  Chapters 
III  and  IV.  However c  the  procedures  presented  here  will  allow  future 
construction  of  simple,  physically  more  consistent,  models  of  interaction 
between  evaporation  and  the  atmosphere. 

2.  Basic  development 

We  begin  with  the  bulk  aerodynamic  relationships  for  surface  moisture  and 
temperature  flux: 

(w'q')  .  -  C  u(q  ,  -  q)  (1) 

'sfc  q  ^sfc 

“  CK  u(T  .  -  T)  (2) 


where  u,  q,  and  T  are  the  atmospheric  wind  speed,  specific  humidity  and 


temperature,  respectively,  measured  at  a  standard  level  such  as  2m.  Cq  and 
Ch  are  non-dimensional  exchange  coefficients,  the  subscript  "sfc"  refers  to 
surface  values,  w  is  the  vertical  motion  and  primes  indicate  turbulent 
fluctuations.  Relationships  (1)  and  (2),  in  principle,  assume  that  the  mean 
flow  is  sufficiently  horizontally  homogeneous  so  that  the  turbulence  is  uniquely 
in  equilibrium  with  the  mean  flow. 

The  exchange  coefficient  appearing  in  the  bulk  aerodynamic  relationship  can 
vary  by  several  factors  with  only  modest  diurnal  variations  of  atmospheric 
stability.  The  same  can  be  said  of  the  coefficient  appearing  in  Dalton's  law 
used  to  derive  the  original  Penman  formula.  In  fact,  the  bulk  aerodynamic 
relationship  and  Dalton's  law  are  essentially  equivalent  (Appendix). 

The  potential  evaporation  can  be  defined  by  replacing  qgfC  in  (1)  with 
the  saturation  surface  value  q*a fc  corresponding  to  the  temperature  of  the 
surface,  in  which  case 

(w'q’)8fc  *  cq  u(q%fc  -  <0-  <3> 

In  analogy  with  the  usual  Penman  derivation,  we  expand  (3)  so  that 

(Wqt)sfC  *  Cq  “  q*)  +  (q*  '  q)1  {4) 

where  q*  is  the  saturation  value  of  the  atmospheric  specific  humidity  at  the 
standard  observational  level. 

To  continue  the  analogy,  we  express  the  saturation  specific  humidity  as  a 
function  of  temperature.  Using  the  relationship 

q*  *  . 


622e*/P 


we  obtain  the  approximation 


.622  de* (T) 


(Tsfc  ‘  T) 


where  e*  is  the  saturation  vapor  pressure  and  p  is  atmospheric  pressure. 
Substituting  (5)  into  (4),  we  obtain 


("^T,s£c  '  Cq  “I 


^  <’,fc  ‘  ”  *«*-«]•  <« 


Using  (2),  the  surface  energy  balance  can  be  written  as 


-p  L  w'q*  +  R  -  S  -  pc  C.  u(T  ,  -  T)  *  0 

v  ^  n  p  h  sfc 


where  p  is  the  surface  density,  Ly  the  specific  latent  heat,  cp  the  specific 
heat  capacity,  R„  the  net  radiative  energy  gained  by  the  surface  and  s  the 
flux  of  heat  to  the  soil  or  vegetation.  Rj,  and  S  are  expressed  in  watts/m2. 
To  eliminate  surface  temperature  we  combine  (6)  and  (7)  and  obtain 


.  (R_  -  S)  pi  C  u(q*  -  q) 


Vch)A) 


PL  (w'q1)  , 
v  sfc 


.622  v  de* (T) 


and  E  is  the  potential  latent  heat  flux  at  the  surface  in  watts/m2. 


Equivalently,  the  evaporation  of  water  in  nm/day  is  3.46xl0“2  E.  The  surface 

-  !  1o3g  A 

moisture  flux  w'q'  in  ms-i  g/kg  is  *  3.1x10  *  E. 

pi» 

V 

With  no  other  information,  we  assume  that  the  turbulent  exchange 
coefficients  for  heat  and  moisture  are  equal.  Then  (8)  becomes 


A(R  -  S)  +  pL  C  u(q*  -  q) 
n  v  g 


1  +  A 


1  +  A 


Relationship  (9)  is  similar  to  other  Penman  formulations  except  that  the  usual 
wind  function  f(u)  is  replaced  by  CqU. 

The  second  term  is  normally  referred  to  as  the  advection  term  since  with  no 
mean  wind  speed  and  no  evaporation  the  specific  humidity  may  approach  satura¬ 
tion.  However,  even  in  the  theoretical  limit  of  vanishing  wind  speed,  turbu¬ 
lence  generated  by  any  surface  heating  will  mix  drier  air  downward,  keeping  air 
near  the  surface  unsaturated.  As  in  the  original  Penman  relationship,  the 
second  term  in  (9)  does  not  vanish  with  vanishing  wind  speed  if  the  air  is  not 
saturated  and  if  the  dependence  of  Cq  on  stability  is  appropriately  chosen. 

The  same  can  be  said  of  the  second  term  when  it  is  identified  with  the  evapora¬ 
tion  measured  by  a  Piche  atmometer  (Brochet  and  Gerbier,  1972).  Even  as  the 
wind  speed  vanishes,  convectively  driven  turbulence  can  ventilate  the  atmometer 
leading  to  evaporation.  For  lack  of  a  better  term,  we  will  refer  to  the  second 
term  as  the  "aerodynamic"  term  although  it  must  be  remembered  that  both  terms  in 
(9)  originate  from  the  bulk  aerodynamic  relationship. 


The  original  Penman  relationship  is  derived  in  the  same  manner  as  (9) 


except  Dalton's  law  (Appendix)  is  used  as  a  starting  point  in  place  of  (1). 
This  relationship  is  typically  expressed  in  the  form  (e.g..  Brutsaert,  1982) 

E  =  4 —  —  +  f(u)  (e*-e)  (10) 

a+y  lv  A+y 


f(u)  =  .26(1+. 54u) 


de*(T) 


c  P 

y  ,  — E — 
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where  E  here  is  expressed  in  un/day,  Q  is  the  net  radiation  flux  density  less 
soil  heat  flux,  and  e*  and  e  are,  respectively,  the  saturation  and  actual  values 
of  vapor  pressure  at  2m.  The  wind  function  f(u)  was  determined  from  evaporation 
pan  measurements  (Penman,  1948).  Many  modifications  of  the  Penman  wind  function 
have  been  suggested  although  the  original  form  still  enjoys  widespread  usage. 
Comparison  of  (10)  and  (9)  indicates  that  the  wind  function  is  proportional  to 
the  exchange  coefficient 

f (u)  «  CqU 

where  the  coefficient  of  proportionality  depends  on  the  units  employed  in  (1) 
and  (10)  (Brutsaert,  1982). 

3.  Dependence  of  exchange  coefficient  on  stability 

The  dependence  of  the  exchange  coefficient  Cq  on  atmospheric  stability 


can  be  expressed  in  terms  of  a  Richardson  number  of  the  form 


Ri  =»  g/0 


~2~ 


u 


where  g  is  the  acceleration  of  gravity,  z  is  the  height  of  the  atmospheric 
observations,  6  is  the  atmospheric  potential  temperature  at  z  and  0sfc  is  the 
potential  temperature  of  the  air  at  the  lower  reference  level.  The  application 
of  Monin-Obukhov  similarity  theory  to  the  bulk  aerodynamic  relationship  requires 
integration  between  two  reference  levels.  The  lower  reference  level  is 
typically  chosen  to  be  the  roughness  height  which  simplifies  the  bulk 
aerodynamic  relationship.  The  derivation  of  the  Penman  relationship  demands 
integration  from  the  surface  where  saturation  is  assumed  in  order  that  the  vapor 
pressure  can  be  determined  from  the  temperature.  The  appropriate  integration 
constant  is  then  the  roughness  length  for  moisture  (Brutsaert,  1982).  Because 
of  the  approximate  nature  of  this  development  and  most  applications  to  actual 
data,  we  will  not  distinguish  between  the  roughness  lengths  for  momentum  and 
scalar  quantities.  For  the  present  data  analysis,  the  influence  of  water  vapor 
on  buoyancy  is  generally  small  and  therefore  also  neglected. 

Based  on  previous  observations  and  certain  asymptotic  constraints,  the 
development  in  Louis  (1979)  along  with  modifications  in  Louis  et  al.  (1982)  lead 
to  the  following  dependence  for  the  unstable  case  (Ri<0) 


where 


(11a) 


75  k 


M1/2 


k  -  .4 
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and  stable  case  (Rl>0) 


(1  +  15Ri) ( 1  +  5Ri)1/2 


(lib) 


Both  formulations  reduce  to  the  usual  log  law  as  Ri+0 .  Note  that  the  exchange 
coefficient  increases  with  increasing  instability  (increasing  negative  Richard¬ 
son  number).  In  the  limit  of  extreme  instability  (Ri-*— »),  after  substituting 
for  the  Richardson  number  into  (11),  the  wind  function  becomes 


C  u  ♦  (8  .  -8)z ] 

q  15  '■z  +  z  '  L8  sfc  J 


Thus,  the  evaporation  rate  becomes  independent  of  wind  speed  and  depends  on 
surface  heating  through  a  square  root  dependence  on  the  surface  temperature 
excess . 

In  the  free  convection  limit,  the  roughness  length  becomes  a  somewhat 
arbitrary  lower  limit  to  the  integration,  which  allows  smooth  matching  to  the 
usual  Monin-Obukhov  similarity  theory.  Ideally  the  expected  rapid  increase  of 
0gfc-9  with  increasing  z/z0  is  such  that  the  free  convection  limit  (12) 
becomes  independent  of  the  roughness  length.  However,  the  actual  sensitivity  of 
(12)  to  the  roughness  length  cannot  be  determined  in  practice  since  the  air 
temperature  becomes  extremely  inhomogeneous  at  the  surface. 

Other  attempts  to  include  the  free  convection  limit  involve  additional 
criteria  and  separate  formulation  of  the  free  convection  case.  These  criteria 
could  be  added  to  the  above  development.  However,  the  free  convection  limit  is 
not  usually  of  practical  importance.  For  example,  in  the  Wangara  experiment. 
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-Ri  rarely  exceeded  one.  As  is  evident  from  Fig.  1  of  Louis  (1979)  for  such 
stability  limits,  the  modelled  influence  of  roughness  length  on  the  stability 
correction  to  the  exchange  coefficient  (11)  is  well-behaved  in  that  it  closely 
approximates  the  original  fit  to  data  by  Businger  et  al .  (1971).  Consequently, 
for  simplicity,  we  proceed  with  (11)  without  further  modification. 

The  stability  corrections  in  the  above  exchange  coefficient  for  heat  and 
water  vapor  is  different  than  that  for  momentum.  However,  the  differences 
between  the  neutral  values  of  the  exchange  coefficients  for  momentum  and  heat  or 
moisture  have  been  neglected  in  the  above  formulation  based  on  Louis  et  al. 
(1982).  This  is  in  contrast  to  the  neutral  limits  in  Louis  (1979)  and  others 
where  the  exchange  coefficient  for  heat  was  larger  than  momentum  and  also 
contrasts  with  Stewart  and  Thom  (1973)  where  the  relationship  between  the 
exchange  coefficients  was  more  complicated. 

4 .  Asymptotic  cases 

We  now  identify  various  special  or  limiting  cases  where  the  wind  speed, 
exchange  coefficient,  humidity  deficit  and/or  net  radiation  less  soil  heat  flux 
vanish  (Table  1).  Combinations  which  do  not  satisfy  the  surface  energy  balance 
have  been  eliminated.  As  before,  molecular  diffusion  of  vapor  is  not 
considered.  The  type  of  evaporation  has  been  classified  as  free  convection  if 
the  wind  speed  vanishes  and  the  implied  turbulence  and  vapor  transport  are 
driven  only  by  buoyancy.  Conversely,  the  evaporation  is  classified  as 
mechanical  if  the  sensible  heat  flux  is  zero  or  downward  in  which  case  vapor  is 
transported  only  by  shear-generated  turbulence. 
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Table  la.  Limiting  cases  of  potential  evaporation  for  the  original  Penman  rela¬ 
tionship. 


Case 

Surface 

Energy 

Balance 

Type 

of 

asaBttiBttv 

Latent  Heat  Flux  E 

Sensible 

Evaporation 

* 

Heat 

Flux 

t  .  ^ 

»  /. 

V- 

■ 

Rad.  term 

aerodynamic 

term 

«  - 

■ 

R  -S 

u 

C 

q*-q 

aero. 

contrib. 

E 

• 

la 

0 

0 

>0 

0 

0 

0 

0 

none 

</• 

2a 

<0 

0 

>0 

0 

0 

<0 

<0 

dewfall 

3a 

>0 

0 

>0 

0 

0 

>0 

>0 

free  conv. 

i 

4a 

>0 

>0 

0 

0 

0 

0 

none 

1 

5a 

>0 

>0 

0 

0 

<0 

<0 

dewfall 

6a 

>0 

>0 

0 

0 

>0 

>0 

:s? 

7a 

0 

>0 

>0 

>0 

>0 

<0 

(see  text) 

8a 

<0 

0 

>0 

>0 

>0 

(1) 

<0 

(see  text) (1) 

’.V 

■.V 

9a 

>0 

0 

>0 

>0 

>0 

>0 

>0 

free  conv. 

10a 

0 

>0 

>0 

>0 

>0 

>0 

<0 

mech.  turb. 

11a 

<0 

>0 

>0 

>0 

>0 

(1) 

<0 

(1) 

12a 

>0 

>0 

>0 

>0 

>0 

>0 

<0 

(2) 

1 

13a 

14a 

>0 

>0 

>0 

>0 

>0 

>0 

>0 

>0 

>0 

>0 

>0 

>0 

0 

>0 

(2) 

(2) 

I  Table 

lb.  Limiting 

cases  for  the  modified  Penman  relationship. 

% 

lb 

0 

1 

1  1 

0 

0 

0 

0 

none 

* 

2b 

<0 

1! 

■9 

0 

0 

<0 

<0 

dewfall 

3b 

>0 

■9 

>0 

0 

0 

>0 

>0 

free  conv. 

fc‘ 

4b 

0 

>0 

>0 

0 

0 

0 

0 

none 

ft. 

5b 

<0 

>0 

>0 

0 

0 

<0 

<0 

dewfall 

6b 

>0 

>0 

>0 

0 

0 

>0 

>0 

ski 

flu 

7b 

0 

D 

■1 

>0 

0 

0 

0 

none 

8b 

<0 

Kfl 

11 

>0 

0 

<0 

<0 

dewfall 

9b 

>0 

19 

>0 

>0 

>0 

>0 

>0 

free  conv. 

£ 

10b 

0 

>0 

>0 

>0 

>0 

>0 

<0 

mech.  turb. 

;X; 

lib 

<0 

>0 

>0 

>0 

>0 

(1) 

<0 

(1) 

12b 

>0 

>0 

>0 

>0 

>0 

>0 

<0 

mech.  turb. 

13b 

>0 

>0 

>0 

>0 

>0 

>0 

0 

mech.  turb. 

m 

14b 

>0 

>0 

>0 

>0 

>0 

>0 

>0 

(1)  depends  on  magnitude  of  each  contribution. 

(2)  12a-13a-14a  because  the  original  relationship  is  independent  of  stability. 


Both  the  original  and  modified  expressions  agree  on  the  existence  of 
evaporation  or  dewfall  for  the  various  cases,  except  for  cases  7  and  8.  Since 
the  original  Penman  wind  function  contains  no  dependence  on  stability,  its  aero¬ 
dynamic  term  incorrectly  predicts  free  convection  of  water  vapor  away  from  the 
surface  for  vanishing  wind  speed  under  stable  conditions  (7a  and  8a)  when  there 
should  be  no  turbulence.  The  same  problem  would  theoretically  occur  with 
vanishing  wind  speed  and  exactly  neutral  conditions.  With  the  stability- 
dependent  exchange  coefficient,  free  convection  of  water  vapor  correctly  occurs 
only  for  unstable  conditions  (upward  heat  flux) . 

5.  Diurnal  variations  in  Wangara 

We  now  compute  the  diurnal  variations  of  potential  evaporation  from  micro- 
meteorological  data  collected  during  the  Wangara  experiment  near  Hay,  Australia 
in  the  winter  of  1967  (Clarke  et  al.,  1971).  The  diurnal  variation  of  stability 
during  the  Wangara  experiment  was,  on  the  average,  less  than  most  would  expect. 
Except  for  day  33,  the  magnitude  of  the  afternoon  Obukhov  length  was  generally 
greater  than  10m  and  on  many  afternoons  greater  than  100m.  This  only  modest 
instability  is  due  to  relatively  low  winter  sun  angles  and  generally  significant 
airflow. 

Potential  evaporation  is  calculated  from  Wangara  data  using  both  the  origi¬ 
nal  Penman  equation  and  the  Penman  equation  modified  to  include  the  stability- 
dependent  exchange  coefficient.  The  40  days  of  data  provided  by  the  Wangara 
program  allow  nearly  900  hourly  calculations  of  potential  evaporation.  Unfor¬ 
tunately,  the  temperature  and  specific  humidity  at  the  reference  height  of  two 
meters  needed  for  the  Penman  calculation  were  not  measured.  These  variables 
were  approximated  by  temperature  and  humidity  data  available  at  a  height  of 
approximately  1.2m.  In  the  daytime,  use  of  the  1.2m  temperatures  would 


overestimate  the  saturation  vapor  pressure  at  2m  while  the  1.2m  specific 
humidity  would  underestimate  the  2m  specific  humidity.  Therefore,  the  net  error 
in  the  estimated  2m  humidity  deficit  is  smaller  than  the  errors  in  the  estimated 
2m  temperature  and  specific  humidity. 

Since  the  surface  temperature  was  not  measured  and  cannot  be  simply  defined 
over  land  surfaces,  the  surface-based  bulk  Richardson  number  could  not  be 
computed.  Instead,  the  layer  Richardson  number  is  computed  using  observations 
at  the  one  and  four  meter  levels.  Because  the  exchange  coefficient  is  a  slowly 
varying  function  of  the  Richardson  number,  except  near  neutral  stability,  the 
error  in  the  estimation  of  the  surface-based  bulk  Richardson  number  will  normal¬ 
ly  cause  much  smaller  errors  in  the  exchange  coefficient.  Note  that  we  cannot 
reintegrate  the  similarity  theory  between  lm  and  4m  to  obtain  a  new  exchange 
coefficient  relationship  because  the  Penman  relationship  demands  integration 
from  the  surface;  that  is,  the  bulk  aerodynamic  relationship  for  moisture  must 
be  defined  with  respect  to  surface  properties  so  that  saturation  can  be  assumed 
allowing  surface  vapor  pressure  to  be  related  to  surface  temperature.  He  assess 
the  importance  of  these  potential  errors  in  Section  8  where  the  exchange  coeffi¬ 
cients  are  also  computed  using  an  iterative  procedure. 

To  compute  the  typical  diurnal  variations,  parameters  for  each  hour  are 
averaged  over  all  of  the  Hangara  days  including  both  sunny  and  cloudy  days. 

Since  p  and  Ly  vary  by  only  a  small  percentage  during  the  day,  they  are 
considered  constant  and  set  equal  to  1.275  kg/m3  and  2.5x10*  J/kg,  respec¬ 
tively.  The  roughness  length  is  assigned  to  be  1.2mm  (Clarke  et  al. ,  1971). 
Figs.  1-2  show  the  diurnal  variation  of  the  remaining  variables  averaged  over  40 
days.  Occasional  missing  observations  contribute  to  some  of  the  hour-to-hour 
noise.  For  most  hours,  less  than  10%  of  the  observations  of  a  given  variable 
were  missing.  As  expected,  the  stability-dependent  exchange  coefficient 
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Fig.  1.  (b)  Mean  diurnal  variation  of  wind  speed  (dot-dash),  specific  humidity 
deficit  (broken),  Louis  stability-dependent  exchange  coefficient 
Cq  (solid),  and  temperature-dependent  coefficient  1/(1+A) 

(dotted) . 
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Averaged  diurnal  variation  of  the  exchange  coefficient  as  computed 
from  the  Louis  formulation  (thick  dot-dash),  and  as  inferred  from 
the  original  Penman  wind  function  f(u)/u  (thick  broken).  Also 
shown  are  averages  for  the  nine  days  with  significant  afternoon 
instability  (thin  dot-dash  and  thin  broken). 


increases  in  the  morning  to  a  maximum  value  occurring  in  the  early  afternoon, 
dropping  off  rapidly  later  in  the  afternoon  to  an  almost  time- independent 
nocturnal  value.  The  inferred  Penman  exchange  coefficient  Cqp  =  f(u)/u,  where 
f(u)  is  the  original  Penman  wind  function,  reaches  a  minimum  in  the  afternoon 
violating  physical  expectations. 

Values  are  also  averaged  for  days  on  which  significant  instability  (Obukhov 
length  <10m)  occurred  in  the  afternoon.  Nine  such  days  are  found.  On  these 
days,  the  afternoon  exchange  coefficient  exceeds  that  inferred  from  the  Penman 
relationship  by  almost  a  factor  of  two  (Fig.  2). 

The  diurnal  variation  of  wind  function  f(u)*CqU  corresponding  to  the 
stability-dependent  aerodynamic  expression  exhibits  significantly  greater 
diurnal  variation  them  the  wind  function  of  the  unmodified  Penman  relationship 
even  when  averaged  over  all  days  (Fig.  3).  Here  the  Penman  wind  function 
follows  a  diurnal  pattern  close  to  that  corresponding  to  that  wind  function  with 
a  constant  neutral  value  of  the  exchange  coefficient,  but  with  a  smaller 
decrease  at  night. 

Fig.  4  shows  the  diurnal  variation  of  the  radiation  term  and  the  various 
aerodynamic  terms.  The  radiation  expression  peaks  around  noon,  whereas  the 
aerodynamic  expressions  peak  in  early  afternoon.  The  aerodynamic  terms  are  as 
large  or  nearly  as  large  as  the  radiation  term  in  contrast  to  some  unstable 
summertime  cases  where  the  radiation  term  is  significantly  larger  than  the  aero¬ 
dynamic  term. 

As  expected,  the  aerodynamic  term  using  the  stability-dependent  exchange 
coefficient  Cq  exhibits,  on  the  average,  considerably  more  diurnal  variation 
than  the  aerodynamic  term  of  the  original  Penman  relationship.  The  aerodynamic 
term  of  the  original  Penman  relationship  is  similar  to  that  aerodynamic  term 
with  a  constant  neutral  exchange  coefficient.  Averaged  over  all  days,  the 
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The  averaged  diurnal  variation  of  the  wind  function  as  computed 
from  the  Louis  stability-dependent  exchange  coefficient  (dot-dash) 
the  original  Penman  formulation  (broken)  and  constant  exchange 
coefficient  (dotted). 
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The  averaged  diurnal  variation  of  the  aerodynamic  term  as  computed 
from  the  Louis  stability-dependent  exchange  coefficient  (dot-dash) 
the  unmodified  Penman  aerodynamic  term  (broken)  and  the  radiation 
term  (dotted) . 


stability-dependent  Penman  relationship  predicts  the  afternoon  aerodynamic  term 
to  be  about  40%  greater  than  that  of  the  original  Penman  relationship  in  the 
afternoon  and  about  50%  less  than  that  of  the  original  Penman  relationship 
during  the  night. 

6.  Nonlinear  diurnal  dependence 

In  most  applications  of  the  Penman  relationship,  the  daily  total  evapora¬ 
tion  is  estimated  by  using  daily  averages  of  wind  speed,  net  radiation,  humidity 
deficit,  and  temperature  and  neglecting  the  influence  of  atmospheric  stability. 
Such  calculations  will  incur  errors  due  to  neglect  of  stability  influences  and 
due  to  the  nonlinear  interaction  (correlation)  between  the  diurnal  variation  of 
variables  which  appear  as  products  in  the  Penman  relationship.  For  example, 
Jobson  (1972)  found  that  nonlinear  interaction  between  the  diurnal  variation  of 
vapor  pressure  and  wind  speed  could  occasionally  lead  to  significant  errors  with 
use  of  da ily-ave raged  variables  in  Dalton's  relationship.  However,  on  90%  of 
the  days,  such  an  error  was  found  to  be  less  than  10%.  This  error  remains 
generally  less  than  25%  even  when  monthly-averaged  variables  are  used  (Hage, 
1975). 

The  results  of  Doorenbos  and  Pruitt  (1977),  also  discussed  by  Stigter 
(1980),  indicate  that  use  of  the  Penman  relationship  for  predicting  24-hour 
evapotranspiration  from  a  well-watered  grass  reference  crop  could  incur  errors 
of  50%  or  more.  The  Penman  relationship  was  found  to  usually  underestimate 
water  loss  with  strong  wind  speed  and  weak  radiative  heating.  This  underestima¬ 
tion  decreased  or  reversed  to  an  overestimation  with  weak  wind  speeds  and  strong 
surface  radiative  heating.  Although  stomatal  resistance  may  have  been  a  factor, 
the  above  variation  of  errors  is  consistent  with  the  influence  of  atmospheric 


stability  on  potential  evaporation  through  the  stability-dependent  exchange 
coefficient  and  its  correlations  with  other  variables  in  the  aerodynamic  term. 

We  now  consider  the  difference  between  calculating  daily  potential 
evaporation  from  daily-averaged  variables  versus  summing  potential  evaporation 
values  computed  from  hourly  variables.  We  refer  to  the  first  case  as 
"linearized  potential  evaporation"  and  the  latter  as  "integrated  potential 
evaporation."  The  original  Penman  relationship  without  stability  influences  is 
relationship  (10)  after  converting  from  vapor  pressure  to  specific  humidity 
using  relationship  (A-3)  in  the  Appendix. 

Table  2  summarizes  the  differences  between  different  expressions  averaged 
over  the  40  days  of  the  Wangara  experiment.  We  first  note  that  the  linearized 
radiation  term  averages  13%  less  than  the  integrated  radiation  term.  This 
underestimation  is  due  to  correlation  between  the  diurnal  variation  of  the 
radiation  (less  soil  heat  flux)  and  the  temperature-dependent  coefficient  of  the 
radiation  term. 

Table  2.  The  daily  aerodynamic  term  averaged  over  the  40  Wangara  days  (mm). 


Method 

Averaged  Value 

%  Difference 
from  Integrated 
Modified 

Averaged 

Absolute 

Difference 

Linearized  Penman 

1.35 

-4% 

15% 

Integrated  Penman 

1.42 

+1% 

12% 

Linearized  Stability- 

Dependent 

0.89 

-36% 

36% 

Integrated  Stability- 

Dependent 

1.40 

— 

— 
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The  linearized  and  daily-integrated  versions  o£  the  aerodynamic  terms  <>t 
the  original  Penman  relationship  average  within  5%  of  each  other.  The  average 
absolute  difference  between  the  two  versions  is  only  6%.  Even  the  maximum 
difference  of  the  40  individual  days  is  only  11%.  This  indicates  that  the 
linearization  due  to  use  of  daily-averaged  values  does  not  lead  to  significant 
errors  in  the  original  Penman  relationship. 

However,  the  linearization  of  the  stability-dependent  aerodynamic  term  does 
create  significant  errors.  Here  the  linearized  aerodynamic  term  averages  less 
than  60%  of  the  integrated  aerodynamic  term.  Differences  on  some  individual 
days  with  large  variations  in  stability  exceeded  90%  of  the  integrated  term.  It 
is  concluded  that  use  of  the  daily-averaged  exchange  coefficient  causes  large 
underestimation  of  the  daily  potential  evaporation  as  is  explicitly  shown  in  the 
next  section.  We  also  found  that  use  of  a  netural  value  of  the  exchange  coeffi¬ 
cient,  corresponding  to  a  logar ithmetic  wind  profile,  underestimates  the  daily 
evaporation  as  previously  concluded  by  Strieker  and  Brutsaert  (1978). 

The  aerodynamic  term  of  the  integrated  original  Penman  averages,  perhaps 
fortuitously,  lies  within  1%  of  that  of  the  integrated  stability-dependent 
Penman  with  an  absolute  difference  averaging  only  12%.  A  maximum  individual 
daily  difference  of  only  24%  of  the  integrated  modified  term  is  found  on  a  day 
with  a  large  diurnal  variation  in  stability. 

The  daily  values  of  the  original  Penman  relationship  agree  rather  closely 
with  the  more  complete  version  partly  due  to  cancellation  of  underestimation  in 
the  daytime  and  overestimation  at  night.  This  agreement  can  also  be  attributed 
to  the  fact  that  the  estimated  roughness  of  1.2mm  is  close  to  the  1.37imn  thought 
to  be  representative  of  pan  conditions  used  to  calibrate  the  original  Penman 
relationship  (Thom  and  Oliver,  1977).  However,  similarity  of  roughness  lengths 
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may  be  a  lesser  factor  since  the  potential  evaporation  does  not  seem  to  be 
especially  sensitive  to  the  roughness  length  (Thom  and  Oliver,  1977). 

Similarly,  the  linearized  aerodynamic  term  of  the  original  Penman  equation 
approximates  the  aerodynamic  term  of  the  integrated,  stability-dependent  Penman 
relationship  reasonably  well.  The  difference  between  them  averages  only  4%  and 
the  absolute  difference  averages  only  15%.  The  maximum  individual  difference, 
about  28%  of  the  integrated  term,  also  occurs  on  a  day  with  a  large  diurnal 
variation  in  stability.  We  conclude  that  for  the  data  considered  here,  the 
original  Penman  relationship  has  been  effectively  calibrated  to  predict  the 
daily  total  evaporation  even  though  it  performs  poorly  on  an  hourly  basis.  We 
further  conclude  that  when  only  24-hour  averages  are  available,  the  original 
Penman  is  preferable  to  the  stability-dependent  Penman  at  least  without  further 
calibration. 

7.  Interactive  terms 

To  study  the  source  of  errors  due  to  use  of  daily-averaged  variables,  each 
variable  has  been  partitioned  into  a  daily  mean  denoted  with  an  overbar  and  an 
hourly  deviation  from  the  daily  mean  denoted  with  a  prime.  Diurnal  variations 
of  density  and  specific  latent  heat  are  much  smaller  than  diurnal  variations  of 
other  parameters  and  are  therefore  assigned  to  be  constant  with  values  of 
p  ■  1.275  kg/m^  and  1^  ■  2.5  x  106  J/kg*  Considering  the  coefficients 
1/(A+1)  and  A/( A+l)  each  to  be  a  single  variable  and  substituting  the  parti¬ 
tioned  variables  into  the  Penman  relationship,  we  obtain 

(1)  (2) 

E  -  +  pLy  (~J)(q*-q)  (a+bu)  +  b(^f)u'  "(q*-q)"'"  + 

(3)  (4)  (5) 

b(q*-q)(^j)T^  +  (a+bu)  (q*-q~  +  bt-j^iVu’  (q*-q)’ 


(13) 
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for  the  original  Penman  relationship  where  a*3.9xlO“3  and  b*2.1xlO”3;  for 
the  stability-dependent  Penman 


(1)  (2) 

\  “  \  +  ®  V<5S=i)  +  (q*"q>’  + 

(3)  (4) 

(X|;rwr  +(^)<«i=5> 

(5)  (6) 


■a+tK'  u’  (q*'q), 

+  cqufe)  <**-*>’  + 

(7) 

(8) 

u* 

*  yirr) '  + 

(9) 

(10) 

u  (q*-q)(”)  Cq* 

+  cq'«i*-q)'  + 

(ID 

n — n - 

(12) 

(q*-q)  (^r) '  Cq'u'  +  (~)  '  Cq’  u'  (q*-qr 


where 

(1)  (2) 

h  ■  Sr  ‘V51  *  Sr)"  ’Vs1' !  • 


(14) 


(15) 


Table  3  summarizes  the  magnitude  of  the  various  terms  averaged  over  the  40 
Wangara  days. 

The  nonlinear  interaction  term  in  the  radiation  expression  averages  about 
13%  of  the  linear  term  averaged  over  all  40  days  and  may  exceed  25%  of  the 
linear  term  on  days  with  large  diurnal  variation  of  temperature  and  net 
radiation  (less  soil  heat  flux). 


Table  3.  Daily  summer  nonlinear  terms  averaged  over  the  40  Wangara  days  (mm) 


Original  Aerodynamic  Term 


averaged-value 

1.45 

0.13 

-0.02 

-0.11 


ratio  to  linear  term 
1.00 
.09 
-.02 
-.08 


Modified  Aerodynamic  Term 


,4 


c 
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In  the  original  Penman  aerodynamic  expression,  the  nonlinear  interaction 
between  wind  and  specific  humidity  deficit  (term  2)  and  between  specific 
humidity  deficit  and  the  temperature-dependent  coefficient  (term  4)  are  found  to 
be  the  most  significant  of  the  nonlinear  terms,  approximately  9%  and  -8%  of  the 
linear  term,  respectively.  On  a  day  with  large  diurnal  variation  of  atmospheric 
variables,  these  terms  reach  about  20%  of  the  magnitude  of  the  linear  term. 

These  two  nonlinear  terms  result  from  higher  wind  speed,  temperature  and 
specific  humidity  deficit  during  the  afternoon  compared  to  nocturnal  periods. 
Note  that  these  two  nonlinear  terms  are  of  opposite  sign  and  approximately 
cancel.  This  explains  why  use  of  the  daily-averaged  values  in  the  original 
Penman  did  not  cause  significant  errors,  at  least  with  data  analyzed  here. 

Of  the  eleven  nonlinear  terms  in  the  stability-dependent  aerodynamic  term, 
seven  are  found  to  be  relatively  unimportant,  (terms  5,  7-12),  summing  to  less 
than  -9%  of  the  linear  term.  The  correlation  between  the  exchange  coefficient 
(Cq)  and  specific  humidity  deficit  (q*-q) ,  and  the  correlation  between  the 
exchange  coefficient  (Cq)  and  wind  (u)  lead  to  the  most  important  nonlinear 
terms  which  average  29%  and  19%  of  the  linear  term,  respectively.  On  a  day  with 
particularly  high  diurnal  variation  of  atmospheric  variables,  these  two 
nonlinear  terms  are  found  to  be  57%  and  26%  of  the  linear  term,  respectively. 
This  strong  correlation  between  diurnal  variations  of  the  stability-dependent 
exchange  coefficient,  wind,  and  humidity  accounts  for  the  large  errors  resulting 
from  the  use  of  daily-averaged  variables  in  the  stability-dependent  Penman 
relationship. 

8.  Iterative  results 

The  Louis  formulation  is  expected  to  incur  errors  associated  with  the 
approximation  of  the  original  similarity  theory.  Errors  also  result  from  the 


use  of  the  layer  Richardson  number  between  one  and  four  meters  in  lieu  of  the 
Richardson  number  evaluated  between  the  surface  roughness  height  and  the 
standard  level  of  2m.  Note  that  in  modelling  situations,  surface  temperature  is 
usually  determined  from  the  surface  energy  balance  and  is  often  quite  different 
from  the  air  temperature  even  if  measured  down  at  the  roughness  height.  This  is 
an  unavoidable  inconsistency  in  modelling  situations  whenever  a  bulk  aerodynamic 
relationship  is  used  in  conjunction  with  a  surface  energy  budget.  That  is, 
similarity  theories  do  not  apply  to  the  actual  surface  temperature  of  the  ground 
even  when  such  temperatures  can  be  defined. 

As  a  check  on  the  modelled  stability  influence  used  here,  we  have  compared 
it  with  the  original  similarity  expressions  of  Businger  et  al.  (1971).  We  have 
integrated  the  similarity  theory  between  the  surface  roughness  height  and  two 
meters  for  wind  and  one  and  four  meters  for  temperature  and  then  used  several 
iterative  approaches  cited  in  the  Introduction.  The  Louis  relationships  used 
here  seem  to  lead  to  a  small  over -estimation  of  the  exchange  coefficient  in 
unstable  situations  during  the  Wangara  experiment,  although  this  disagreement 
varied  somewhat  with  the  choice  of  iterative  scheme.  In  addition,  the  original 
similarity  expressions  are  uncertain  in  cases  of  large  instability  or  large 
stability.  Lower  limits  on  the  Obukhov  length  must  be  imposed  for  stable  situa¬ 
tions  in  order  to  insure  physically  realistic  behavior  and  in  some  cases  insure 
convergence  of  the  iterative  scheme. 

We  conclude  that  the  qualitative  differences  between  the  original  Penman 
relationship  and  those  modified  to  include  stability  dependencies  are  not 
critically  dependent  on  the  stability  formulation  as  also  concluded  by  Strieker 
and  Brutsaert  (1978).  However,  even  the  modified  formulations  are  only  approxi¬ 
mate  and  a  precise  quantitative  evaluation  of  the  original  Penman  relationship 


9.  Radiation -Richard son  number 


The  iterative  procedure  may  be  too  cumbersome  for  many  applications  while 
the  evaluation  of  the  Richardson  number  requires  temperature  at  two  levels,  not 
typically  available  in  modelling  and  routine  observational  situations.  The 
Richardson  number  can  be  modified  by  replacing  the  temperature  difference  with  a 
dependence  on  radiation.  The  resulting  parameter  is  more  "external"  than  the 
usual  Richardson  number,  since  the  temperature  difference  and  turbulent  fluxes 
are  directly  coupled.  The  surface  heat  flux  is  uniquely  related  to  net  radia¬ 
tion  under  conditions  of  potential  evaporation  and  negligible  heat  flux  to  the 
soil.  However,  in  general,  the  evaporation  is  less  than  potential  and  the  ratio 
of  the  heat  flux  to  the  surface  net  radiation  varies. 

A  "radiation-Richardson"  number  can  be  developed  by  beginning  with  the  flux 
Richardson  number 

Rf  =  (16) 

jbr  w'y* 

3z 


where  w'61  is  the  surface  temperature  flux,  w'v'  is  the  surface  momentum  flux 
and  v  is  the  mean  flow  vector.  He  replace  the  seldom-measured  surface  heat  flux 
with  the  radiative  temperature  flux  less  soil  heat  flux 


R 


(R  -S) 
n 


(17) 


Rn  is  positive  with  net  downward  radiative  flux  and  is  expressed  in  watts/m2 
in  which  case  R  is  in  units  of  K  ms-1.  In  most  practical  situations,  S  would 
be  neglected.  Scaling  velocity  fluctuations  with  u,  the  wind  speed  at  height  z, 
a  scale  value  for  the  flux  Richardson  number  (16)  becomes  proportional  to  the 
radiation-Richardson  number 

RiRad  =  (g/0)  R  z/u3  •  (18) 

The  relationship  between  the  radiation-Richardson  number  and  atmospheric 
stability  depends  on  the  actual  evaporation  rate  so  that  the  radiation-Richard¬ 
son  number  is  only  a  crude  estimate  of  atmospheric  stability.  However,  since 
the  main  variation  of  the  exchange  coefficient  occurs  in  the  transition  between 
stable  and  unstable  cases,  the  crude  estimate  of  stability  based  on  the  radia¬ 
tion-Richardson  number  will  be  of  utility. 

To  develop  the  intended  use  of  the  radiation-Richardson  number,  we  regress 
z/L  on  the  radiation-Richardson  number  where  L,  the  Obukhov  length,  is  computed 
iteratively  using  the  similarity  relationships  of  Businger  et  al.  (1971).  The 
height  z  is  again  2m.  Cases  where  net  radiative  heat  gain  is  positive  and  z/L 
is  stable  and  vice  versa  are  eliminated  since  these  cases  normally  occur  in 
transitional  periods  when  similarity  theory  is  not  expected  to  apply.  Fortu¬ 
nately,  potential  evaporation  rates  are  small  during  such  periods. 

The  radiation-Richardson  number  and  z/L  are  linearly  correlated  with  a 
coefficient  of  .95  for  the  unstable  cases  and  .57  for  stable  cases.  However, 
distributions  of  these  parameters  are  strongly  skewed.  The  cube  root  of  z/L  and 
the  radiation-Richardson  number  are  more  normally  distributed  and  will  be  used 
for  the  regression  relationships.  Note  that  (z/L) 1/3  is  inversely  related  to 
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the  suiface  friction  velocity  while  (Ripft<<)  *-/3  is  inversely  proportional 

to  the  wind  speed,  so  that  the  resulting  regression  relationship  is  analogous  to 

a  resistance  law. 

The  cube  root  of  z/L  is  correlated  to  (RiRad)3"^3  with  a  correlation 
coefficient  of  .90  in  the  unstable  case  and  .77  in  the  stable  case.  The  regres¬ 
sion  relationships  for  the  stable  and  unstable  cases  are: 


Vi- 


(z/L)V3  »  -  8.64  (RiRad)1/3-^;  stable  case 

(19) 

(z/L) l/3  =  -15.29  (RiRad)1/3"*13*  unstable  case. 

Both  relationships  predict  that  z/L  approaches  about  -10~3  as  the  net  radia¬ 
tion  vanishes.  This  small  constant  has  no  special  significance  for  the  near 
neutral  case  but  rather  improves  the  fit  over  the  range  of  the  values  of  the 
radiation-Richardson  number.  A  higher  order  model  is  not  justified  because  of 
the  very  approximate  nature  of  this  development. 

The  above  regression  relationships  (19)  were  used  to  estimate  z/L  and 
subsequently  to  compute  the  exchange  coefficients  for  the  Wangara  data.  These 
exchange  coefficients  averaged  over  the  forty  days  appear  to  lead  to  an  under¬ 
estimation  of  the  stability  influence  (Fig.  5).  With  wetter  surface  conditions, 
this  technique  may  overestimate  the  stability  influence.  However,  these  simple 
explicit  relationships  based  on  the  radiation-Richardson  number  should  be  a 
significant  improvement  upon  complete  neglect  of  the  influence  of  atmospheric 
stability  and  at  the  same  time  do  not  require  additional  observations  as  with 
other  stability  parameters. 


10.  Penman-Monteith 

To  include  the  influence  of  stomatal  control,  the  Penman  relationship  is 
often  multiplied  by  a  plant  coefficient  which  is  generally  less  than  unity.  As 
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TIME  OF  DAY  (hrs) 

The  averaged  diurnal  variation  of  the  exchange  coefficient  based  on 
the  Louis  formulation  (dot -dash) ,  the  radiation-Richardson 
number-regression  relationship  (solid)  and  as  inferred  from  the 
Penman  wind  function  (broken).  Also  shown  is  the  atmospheric 
resistance  coefficient  (thin-broken). 


an  alternative,  the  Penman-Monteith  relationship  (Monteith,  1965)  is  frequently 
used.  This  approach  introduces  the  influence  of  vegetation  properties  through  a 
surface  resistance  factor  rs  and  the  aerodynamic  exchange  coefficient  is 
absorbed  into  a  coefficient  for  aerodynamic  resistance  to  atmospheric  vapor  flux 
ra  through  the  relationship 

ra  =  1/ ( Cqu )  (20) 

where  Cq  is  parameterized  according  to  (11).  Thus  the  resistance  coefficient 
is  sensitive  to  atmospheric  stability.  We  have  computed  the  diurnal  variation 
of  the  resistance  coefficient  from  the  exchange  coefficient  averaged  over  all 
the  days  for  each  hour.  The  resistance  coefficient  cannot  be  averaged  directly 
since  during  very  stable  conditions,  the  resistance  coefficient  becomes  orders 
of  magnitude  greater  than  typical  values  and  theoretically  can  become  infinite. 

Diurnal  variations  of  ra,  as  computed  from  the  Wangara  data,  are  plotted 
in  Fig.  5.  Since  the  diurnal  variation  of  ra  is  substantial ,  the  neglect  of 
stability  and  nonlinear  interactions  in  the  Penman-Monteith  expression  is 
expected  to  lead  to  large  errors.  The  diurnal  variation  of  ra  could  in  turn 
be  used  to  assess  the  diurnal  variation  of  the  coefficient  a  in  the  Priestly- 
Taylor  model  by  employing  the  model  of  de  Bruin  (1983). 

Relationship  (20)  along  with  either  (11),  (19)  or  the  procedure  in  Section 
8  would  allow  inclusion  of  the  stability  influence  in  the  Penman-Montieth 
expression. 

11.  Conclusions 

The  potential  evaporation,  as  related  exclusively  to  atmospheric  variables, 
is  found  to  be  quite  sensitive  to  the  diurnal  variation  of  the  exchange 
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coefficient  appearing  in  the  bulk  aerodynamic  relationship  (Dalton's  law).  The 
usual  neglect  of  such  diurnal  variations  of  stability  leads  to  a  factor  of  two 
smaller  estimates  of  the  values  for  the  aerodynamic  term  on  afternoons  with 
modest  instabilities  although  the  differences  are  reduced  substantially  when 
averaged  over  all  of  the  days  in  the  Wangara  experiment  (Section  5) .  Strong 
instability  (very  small  Obukhov  length)  did  not  occur  with  the  data  analyzed 
here  due  to  the  relatively  low  winter  sun  angle  and  generally  significant 
airflow.  In  many  summer  mid-latitude  situations,  the  diurnal  variation  of 
stability  will  be  greater  than  found  here.  Over  a  fully  vegetated  surface  with 
moist  soil  and  significant  airflow,  the  diurnal  variation  will  often  be  less 
than  found  here. 

Although  responding  inadequately  to  diurnal  variations,  the  original  Penman 
relationship  predicted  daily  totals  of  the  potential  evaporation  which  are  in 
good  agreement  with  values  predicted  by  the  more  complete  relationship.  This 
can  be  attributed  to  partial  cancellation  of  the  daytime  and  nocturnal 
influences  in  atmospheric  stability  and  the  similarity  between  the  roughness 
length  in  the  Mangara  experiment  and  that  corresponding  to  the  original 
calibration  of  the  Penman  relationship.  Use  of  daily-averaged  variables  in  lieu 
of  summing  hourly  estimates  from  the  original  Penman  relationship  led  to  little 
error. 

We  have  also  developed  a  new  formulation  for  computing  the 
stability-dependent  exchange  coefficient  by  defining  a  "radiation-Richardson 
number”  (Section  9).  Although  less  accurate,  this  simple  formulation  uses  only 
variables  which  are  already  required  for  evaluation  of  the  unmodified  Penman 
relationship.  Variations  of  atmospheric  stability  may  also  lead  to  significant 
day-to-day  and  regional  variations  of  total  potential  evaporation. 


Appendix:  Relationships  between  specific 
humidity  and  vapor  pressure 

The  Penman  relationship  was  originally  derived  and  generally  applied  in 
terms  of  Dalton's  law  which  relates  moisture  flux  to  the  gradient  of  vapor 
pressure.  Turbulent  transport  is  better  related  to  a  conservative  moisture 
variable  such  as  specific  humidity  as  used  in  the  bulk  aerodynamic  approach. 

Here  we  relate  the  two  approaches. 

Specific  humidity  is  related  to  vapor  pressure  through  the  approximate 
relationship 

q  *  ee/p  (Al) 

where  e  is  the  ratio  of  the  molecular  mass  of  vapor  to  that  of  dry  air,  equal  to 
approximately  0.622. 

Changes  in  specific  humidity  are  then 

.  ,M) 

p  p2 

Changes  of  vapor  pressure  between  the  surface  and  some  standard  level,  such  as 
2m,  can  be  O(10~l)e  in  cases  of  significant  moisture  flux.  The  same  can  be 
said  of  fluctuations  of  vapor  pressure  due  to  turbulence.  On  the  other  hand 
turbulent  pressure  fluctuations  near  the  surface  can  be  expected  to  be 
0(10~4)p  so  that  the  ratio  of  the  first  term  to  the  second  term  in  (A2)  is 


large  in  which  case 


The  dependence  of  saturation  specific  humidity  on  temperature  is  related  to 
the  dependence  of  saturation  vapor  pressure  on  temperature  through  (A2)  such 
that 


dq*  £  de»  e*  dp 

dT  p  dT  G  2  dT 
P 


(A4) 


The  ratio  of  the  second  term  to  the  first  one  is 

e£  dp  de* 
p  dT  '  dT 


(A5) 


The  relationship  between  pressure  and  temperature  is  restricted  by  the 
ideal  gas  law  and  thermodynamic  relationship.  In  the  case  of  adiabatic  motion 
of  turbulent  elements  the  relationship  between  temperature  and  pressure  is 


dT 


dp  . 


(A6) 


Using  the  definition  of  potential  temperature 


dT 


Po/Tk 


(A7) 


where  k  ■  .286  and  pQ  -  1000  mb. 

Relationship  (A7)  implies  that  the  ratio  (A5)  is 


dT  .  de* 


(A8) 


Using  the  Claus ius-Clapeyr on  equation. 
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the  ratio  (8)  becomes 


R  T 
v 

icL 


(A10) 


Where  Rv  is  the  specific  gas  constant  for  water  vapor  and  Ly  is  the  latent 
heat  of  condensation.  This  ratio  has  a  maximum  value  of  about  .2.  Then  the 
second  term  in  (A5)  can  be  neglected  for  a  rough  approximation  in  which  case 


dq*  _e  de* 
dT  *  p  dT 


(All) 


Integrating  this  relationship  from  the  surface  temperature  to  the  temperature  at 
some  standard  level,  we  obtain 


i*  —  n* 


sfc 


-  (e*  -  e*  ) 
p  sfc 


(A12) 


where  the  height-dependent  pressure  can  be  replaced  with  the  pressure  measured 
at  the  standard  level  with  an  error  of  only  £10~3.  Integrating  (A3)  at  a 
constant  pressure  level  produces  the  same  relationship  for  the  general  case 
where  saturation  is  not  assumed.  Eq.  (A12)  verifies  that  the  same  approximation 
holds  with  saturation  as  represented  by  the  Clausius-Clapeyron  equation.  Thus 


we  can  write 


which  is  valid  with  or  without  saturation 


In  terns  of  the  exchange  coefficient,  Dalton's  law  can  then  be  expressed  as 


(w'q' ) 


sfc 
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In  the  development  of  this  chapter,  Dalton's  law  is  not  used.  This  appendix 
shorn  that  the  use  of  Dalton's  law  is  equivalent  to  use  of  the  bulk  aerodynamic 
relationship  which  neglects  the  stability  dependence  of  the  exchange  coeffi- 
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III.  A  Two-Layer  Model  of  Soil  Hydrology 
Abstract 

A  two-layer  model  of  soil  hydrology  is  developed  for  applications  where 
only  limited  computer  time  and  complexity  are  allowed.  Volumetric  soil  water  is 
computed  in  a  thin  upper  layer  for  use  in  calculation  of  surface  evaporation. 
Storage  of  water  is  computed  for  an  underlying  deeper  layer. 

In  {in  effort  to  identify  the  influence  of  significant  asymmetric  truncation 
errors  in  the  two-layer  model,  this  model  is  compared  with  the  100-level  model 
of  Boersma  et  al.  (1983).  Comparisons  are  made  for  modelled  soils  with  clay, 
loam  and  sand  properties  for  various  time  dependencies  of  potential  evaporation 
and  precipitation.  Truncation  errors  in  the  resulting  two-layer  model  appear  to 
be  modest  at  least  compared  to  errors  asociated  with  difficulty  in  estimation  of 
the  hydraulic  diffusivity  and  its  strong  dependence  on  soil  water  content. 

Minimization  of  the  influence  of  truncation  errors  requires:  1)  choosing 
the  upper  layer  to  be  sufficiently  thin,  2)  allowing  the  soil  water  gradient  to 
directly  control  surface  evaporation  and  3)  using  suitable  numerical  implementa¬ 
tion  for  the  evaluation  of  internal  soil  water  flux. 
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1.  Introduction 

The  goal  of  this  study  is  to  develop  a  model  of  soil  hydrology  which  is 
simple  enough  to  use  in  concert  with  atmospheric  general  circulation  models, 
forecast  models  and  various  operational  situations  where  computer  time  and 
allowed  complexity  are  restricted.  At  the  same  time  a  certain  minimum  physics 
seems  necessary  to  include  for  most  applications.  For  example,  surface  evapora¬ 
tion  is  related  to  soil  moisture  near  the  surface  while  storage  of  soil  water 
and  transpiration  are  related  to  soil  moisture  in  a  deeper  layer.  Therefore  at 
least  two  layers  including  a  thin  upper  layer  seem  desirable.  This  geometry  was 
proposed  in  Deardocff  (1977)  where  the  thickness  of  the  upper  layer  was  regarded 
as  vanishingly  small. 

The  strong  dependence  of  hydraulic  diffusivity  on  soil  water  content  is  a 
second  necessary  inclusion  for  most  applications.  As  a  result  of  this  depend¬ 
ence,  the  hydraulic  diffusivity  can  vary  by  orders  of  magnitude  with  depth  or 
with  time  at  a  given  level  during  the  diurnal  drying  cycle.  In  the  two-layer 
model  of  Deardorff  (1977),  this  dependence  appeared  to  be  indirectly  included  as 
a  result  of  calibration  of  an  empirical  evaporation  function  using  observations 
of  drying  in  Adelanto  loam  (Jackson,  1973).  This  indirect  calibration  would 
probably  have  to  be  significantly  altered  for  other  soil  types  since  the 
hydraulic  diffusivity  is  also  sensitive  to  soil  type. 

Jersey  (1982)  more  directly  includes  the  dependence  of  hydraulic  diffusiv¬ 
ity  on  soil  water  content  in  a  two-layer  soil  model  by  using  the  basic  soil 
water  transport  relationships.  These  equations  form  the  basis  of  the  high- 
resolution  soil  models  (Nimah  and  Hanks,  1973;  Feddes  et  al.,  1974;  McCumber  and 
Pielke,  1981;  Camillo  et  al. ,  1983;  and  Boersma  et  al. ,  1983). 

A  third  crucial  feature  is  suitable  representation  of  the  relationship 


between  soil  evaporation  and  the  atmospheric  potential  evaporation.  The 


usual  technique  of  relating  evaporation  to  layer  average  soil  water  content  is 
not  very  well  posed.  Use  of  information  on  near-surface  soil  water  flux  is  much 
preferable  since  it  allows  control  by  the  soil  moisture  profile. 

Using  this  modelling  approach  for  evaporation,  we  develop  a  two-layer  model 
of  soil  hydrology  from  the  basic  transport  equations  cited  above.  As  in 
Deardorff  (1977),  the  upper  layer  is  chosen  to  be  thin  so  as  to  better  approxi¬ 
mate  the  large,  near-surface  gradients  associated  with  the  diurnal  drying 
cycle.  The  use  of  basic  transport  equations  in  place  of  the  force-restore 
approach  allows  relating  surface  evaporation  to  basic  characteristics  of  the 
soil. 

The  truncation  errors  associated  with  any  two-layer  model  require  special 
attention  since  they  are  large  and  of  different  nature  from  the  usual  truncation 
errors  in  higher  resolution  models.  Here  we  study  the  behavior  of  truncation 
errors  by  comparing  this  two-layer  model  with  the  high  resolution  model  of 
Boersma  et  al .  (1983).  Hie  depth  of  the  upper  layer  and  the  method  for  imple¬ 
menting  the  numerical  calculation  of  internal  water  flux  significantly  affect 
the  magnitude  of  the  truncation  errors. 

In  this  study,  we  do  not  consider  the  influences  of  vertical  temperature 
gradients  on  soil  water  flux  as  in  Camillo  et  al.  (1983)  and  Jersey  (1982). 
Distinction  is  not  made  between  liquid  transport  and  the  generally  less 
important  vapor  transport.  Since  the  model  is  one -dimensional ,  horizontal 
inhomogeneity  is  excluded.  While  soils  are  virtually  everywhere  inhomogeneous 
in  a  significant  way,  it  is  assumed  that  a  suitable  hydraulic  diffusivity  exists 
for  approximating  vertical  transport  of  water  over  a  reference  area. 


2.  Basic  model 


a.  Layer-integration 

Neglecting  root  uptake,  the  transport  of  water  in  the  soil  can  be  described 
by  the  equation  (e.g.,  Hillel,  1980) 
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where  6  is  the  volumetric  water  content,  D(0)  is  the  soil  water  diffusivity, 

K(  0)  is  referred  to  as  the  hydraulic  conductivity  although  it  has  units  of 
velocity  and  negative  z  is  increasing  downward.  Integrating  (1)  over  the  upper 
layer  of  depth  (Fig.  1)  we  obtain 
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is  the  direct  evaporation  of  soil  water  to  the  atmosphere.  I  is  the 
infiltration  rate  which  is  zero  in  the  absence  of  precipitation  and  otherwise  is 
equal  to  the  precipitation  rate  up  to  a  maximum  infiltration  rate.  This  maximum 
rate  is  parameterized  in  terms  of  the  estimated  flux  at  the  surface  under  condi¬ 
tions  of  saturation,  approximated  as 

•  (3> 

Here  the  gradient  near  the  surface  has  been  estimated  in  terms  of  a  linear 
variation  between  the  saturation  value  at  the  surface  and  a  layer  mean  value 


50 


'V  ■  T -■  .*.■  .y 77  .<•■  \-  .-■  -7V  r»  yr 


assumed  to  occur  at  the  mid-point  o£  the  upper  layer.  Any  precipitation  which 
cannot  infiltrate  or  re-evaporate  is  specified  to  be  runoff.  In  the  case  of 
heavy  precipitation  and  a  large  time  step  such  as  one  hour  or  greater,  it  may  be 
necessary  at  the  end  of  the  time  step  to  remove  any  soil  water  in  excess  of 
saturation  as  additional  runoff. 

We  obtain  for  the  lower  layer  (see  Fig.  1): 
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where  K(0)_Z2  is  the  moisture  loss  due  to  "gravitational"  percolation  through 
the  bottom  of  the  lower  layer.  For  simplicity,  the  diffusion  of  soil  water 
across  the  bottom  of  the  lower  layer  z*-z2  is  neglected  since  gradients  and 
resulting  fluxes  at  this  depth  are  generally  unimportant  on  the  time  scale  of  a 
few  weeks  or  less.  Exceptions  include  a  high  water  table  and  the  advance  of  a 
deep  "wetting  front."  In  the  latter  case,  the  thickness  of  the  model  should  be 
increased. 

b.  Model  geometry 

For  evaluation  of  surface  evaporation,  the  upper  layer  should  be  chosen  as 
thin  as  possible  so  as  to  be  representative  of  the  water  content  at  the 
surface.  However,  assigning  the  layer  to  be  too  thin  has  two  disadvantages. 
First  direct  comparison  with  observations  are  more  difficult  since  vertical 
gradients  and  horizontal  inhomogeneity  may  be  extremely  large  near  the  surface 
making  the  water  content  difficult  to  define  from  observations.  Secondly,  a 
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large  difference  in  thickness  between  the  two  layers  leads  to  asymmetric  finite 
differencing  and  generally  more  complicated  truncation  errors  in  the  approxima¬ 
tion  of  vertical  gradients  between  the  two  layers.  As  a  compromise  between 
these  considerations  and  cognizance  of  usual  observational  levels,  we  tentative¬ 
ly  choose  the  thickness  of  the  upper  layer  to  be  5  cm  and  the  thickness  of  the 
lower  layer  to  be  95  cm. 

Implicit  time  differencing  is  used  with  respect  to  evaluation  of  the 
vertical  gradient.  The  time  step  is  30  minutes. 


c.  Flux  parameterization 

The  hydraulic  conductivity  and  diffusivity  are  not  properties  of  the  soil 
but  rather  strong  functions  of  the  water  content  and  history  of  soil  wetting  and 
drying.  The  rapid  decrease  of  hydraulic  conductivity  and  diffusivity  with 
decreasing  volumetric  soil  water  content  found  in  several  previous  studies  is 
shown  in  Figs.  2-3.  These  figures  show  the  large  differences  between  soil  types 
but  do  not  include  hysteresis  effects  due  to  dependence  on  soil  history. 

We  will  assume  that  the  conductivity  and  diffusivity  behave  in  a 
sufficiently  organized  fashion  such  that  its  dependence  on  water  content  can  be 
approximated  with  simple  functions.  As  in  Boersma  et  al.  (1983)  (hereafter 
referred  to  as  B),  we  parameterize  the  hydraulic  diffusivity  and  conductivity 
with  piecewise  log-linear  fit  to  data  from  Gardner  (1960).  We  have  also 
operated  the  two-layer  model  with  the  hydraulic  relationships  presented  in  Clapp 
and  Hornberger  (1978)  and  Campbell  (1974)  and  applied  in  McCumber  and  Pielke 
(1981).  These  relationships  are  based  on  a  convenient  mathematical  format 
adopted  in  Gardner  et  al.  (1970).  Such  relationships  allow  easy  application  to 
a  variety  of  soil  types.  Here  we  use  the  fit  to  data  from  Gardner  (1960)  to 
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facilitate  comparison  with  B. 
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Examples  of  the  dependence  of  soil  hydraulic  diffusivity  on  volu¬ 
metric  soil  water  content  from  loam,  (HB^,  Hanks  and  Bowers, 
1962);  (J,  Jackson,  1973);  (GHB,  Gardner  et  al . ,  1970);  silt  loam 
(HBg,  Hanks  and  Bowers,  1962);  clay  (P,  Passioura  and  Cowan, 
1968);  results  approximated  from  Gardner  (1960)  for  sand  (Bg), 
loam  (B^)  and  clay  (Be) ;  relationship  from  Clapp  and  Hornberger 
(1978)  for  sand  (CHg) ,  loam  (CH^)  and  clay  (CHc) . 
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VOLUMETRIC  WTER  CONTENT,  9 

Fig.  3.  Examples  of  the  dependence  of  hydraulic  conductivity  on  volumetric 
soil  water  content  for  sand,  (DL,  Day  and  Luthin,  1956);  (Black  et 
al. ,  1970,  0-50  cm-BGTi,  50-150  cm-BGT2) ,  loam  (J,  Jackson, 

1973) ;  (MHLi  and  MHL2,  Marshall  and  Holmes,  1979) ;  (GHB, 

Gardner  et  al. ,  1970) ;  results  approximated  from  Gardner  (1960)  for 
sand  (Bs) ,  loam  (B^)  and  clay  (Be) ;  relationship  from  Clapp 
and  Homberger  (1978)  for  sand  (CHs) ,  loam  (CH^)  and  clay 
(CHC). 
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The  gradient  of  soil  water  at  the  interfacial  level,  -z\ ,  is  approximated 
by  assuming  that  the  layer-averaged  water  content  occurs  at  the  layer  mid-level 
in  which  case  the  present  approach  becomes  equivalent  to  the  usual  finite 
differencing  without  layer-integration.  Assuming  linear  variation  between  the 
two  mid-levels  we  obtain 
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where  Z2/2  is  the  distance  between  mid-levels.  This  approximation  will 
normally  significantly  underestimate  the  gradient  of  soil  water  at  the 
interfacial  level  in  periods  of  evaporation  when  the  gradient  decreases  rapidly 
with  depth. 

Due  to  this  decrease  with  depth,  linear  interpolation  between  mid-levels 
will  also  significantly  underestimate  the  soil  water  content  at  the  interface 
and  therefore  underestimate  the  diffusivity  at  the  interface.  During 
evaporation  periods,  the  change  in  water  flux  with  depth  is  typically  a  small 
difference  between  two  much  larger  effects,  one  being  the  rapid  decrease  of  the 
water  content  gradient  with  depth  and  the  other  being  the  rapid  increase  of 
hydraulic  diffusivity  with  depth.  To  minimize  truncation  errors  due  to  the 
large  distance  between  mid-levels,  it  is  important  that  the  computed  values  of 
the  diffusivity  and  gradient  are  self-consistent  or  represent  approximately  the 
same  level.  This  is  more  important  than  insisting  that  the  estimate  of  the  flux 


coincides  with  the  inter  facial  level.  Because  of  the  curvature  of  the  water 
content  profile,  the  bulk  gradient  (5)  is  representative  of  the  gradient  well 
within  the  lower  layer.  We  found  that  during  drying,  use  of  the  hydraulic 
diffusivity  and  conductivity  corresponding  to  the  water  content  of  the  lower 
layer  instead  of  use  of  interpolated  or  vertically  averaged  values  produced 
results  closer  to  the  high  resolution  model  B  (Section  4) . 

After  periods  of  weak  evaporation  when  the  soil  water  varies  slowly  with 
depth,  the  various  methods  of  estimating  the  interfacial  water  flux  will  produce 
comparable  estimates.  In  the  case  of  precipitation,  the  downward  advance  of  the 
wetting  front  is  determined  more  by  the  hydraulic  properties  of  the  upper  wetted 
soil.  Using  the  diffusivity  and  conductivity  of  the  upper  layer  of  the 
two-layer  model  during  precipitation  events  (Section  5)  appeared  to  minimize 
truncation  errors. 

Therefore,  when  the  lower  layer  is  more  moist,  the  diffusivity  and 
hydraulic  conductivity  at  the  interface  z,  (Fig.  1)  is  evaluated  using  the  soil 
water  content  of  the  lower  layer  02  so  that 

(K  ,  D  )  -  [K ( 0  ) ,  D(0  )];  0,  >  0 ,  .  (6) 

z,  z  2  2  2  i 


When  the  upper  layer  is  more  moist,  the  diffusivity  and  hydraulic 
conductivity  at  the  interface  z  are  evaluated  from  the  water  content  of  the 
overlying  upper  layer  0-j  in  which  case 
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For  the  cases  considered  here,  this  procedure  reduces  truncation  errors  as 
effectively  as  the  finite  element  approach  and  does  not  require  the  more 


complicated  boundary  conditions  needed  for  the  finite  element  approach.  The 
infiltration  condition  (3)  is  consistent  with  (6-7)  in  *nat  the  diffusivity  and 
hydraulic  conductivity  are  evaluated  from  "upstream"  water  content. 

3.  Vertical  structure 

We  first  study  the  bulk  vertical  structure  of  soil  water,  its  dependence  on 
evaporation  and  precipitation  and  resulting  truncation  errors  associated  with 
the  estimate  of  the  interfacial  gradient  (5).  The  magnitude  of  potential 
truncation  errors  in  low  resolution  models  will  be  estimated  using  the  high 
resolution  model  B  where  truncation  errors  are  expected  to  be  small.  In  B,  (1) 
is  evaluated  at  each  of  100  levels.  The  vertical  resolution  varies  from  1/2  cm 
near  the  surface  to  2  cm  near  the  bottom  (158  cm) . 

As  a  single  index  for  potential  truncation  errors  in  low  resolution  models, 
we  define  the  ratio  of  the  local  gradient  at  5  cm  to  the  bulk  gradient  between 
2  1/2  cm  and  52  1/2  cm,  both  computed  from  model  B.  The  local  gradient  at  5  cm 
is  computed  by  differencing  between  the  4  and  6  cm  levels  in  B.  The  ratio  of 
gradients  will  allow  quick  interference  of  the  behavior  of  truncation  errors  due 
to  profile  curvature  of  the  water  content. 

Comparisons  are  made  for  Chino  Clay,  a  soil  with  very  low  hydraulic 
diffusivity  and  Pachappa  sand  which  has  large  hydraulic  diffusivity.  Two  types 
of  experiments  are  run,  one  with  specified  constant  potential  evaporation 
(atmospheric  demand)  of  .1  cm/hour,  typical  of  moderate  daytime  evaporation,  and 
one  where  the  potential  evaporation  is  specified  to  be  constant  (.1  cm/hour)  for 
six  hours,  zero  for  eighteen  hours  and  so  forth.  The  latter  experiments  allow 
examination  of  the  soil  response  to  changes  of  potential  evaporation,  and  when 
viewed  over  several  days,  crudely  simulate  diurnal  variations  of  evaporation. 
Numerical  experiments  were  also  run  for  the  case  of  sinusoidal  variation  of 
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Ratio  of  the  bulk  gradient  of  volumetric  water  content  to  the  local 
gradient  at  5  cm  as  computed  from  the  high-resolution  model.  Ratio 
values  significantly  greater  than  one  reflect  the  failure  of  the 
bulk  gradient  to  recognize  the  strong  curvature  near  the  surface 
which  is  induced  by  surface  evaporation.  Deviations  of  the  ratio 
value  from  unity  indicate  potential  truncation  errors  in  the  two- 
layer  model.  See  text  for  further  explanation. 
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evaporation.  Initial  conditions  are  chosen  to  be  the  same  as  those  in  B;  8*. 29 
for  Pachappa  sand  and  8-. 47  for  Chino  Clay. 

Fig.  4  indicates  that  the  actual  gradient  at  5  cm  is  significantly  larger 
than  the  bulk  gradient  during  evaporation  periods  since  surface  drying  induces  a 
moisture  gradient  which  is  large  at  the  surface  and  decreases  rapidly  with 
depth.  The  ratio  of  the  two  moisture  gradients  decreases  with  time  as  the 
influence  of  surface  drying  "diffuses"  to  lower  levels.  After  about  5  hours 
with  constant  potential  evaporation,  the  ratio  of  the  local  gradient  at  5  cm  to 
the  bulk  gradient  in  the  Pachappa  sand  has  decreased  to  a  nearly  time- 
independent  value  of  about  3.75.  This  ratio  appears  to  approach  the  same 
asymptotic  value  in  the  Chino  Clay  but  requires  a  much  longer  time  due  to  lower 
diffusivity.  Without  considering  errors  in  the  diffusivity,  this  ratio  implies 
that  the  between-layer  flux  computed  in  the  two-layer  model  will  be  underesti¬ 
mated  by  a  factor  of  almost  four. 


In  the  diurnal  simulations,  both  the  daily-averaged  ratio  of  gradients  and 
the  diurnal  amplitude  of  the  ratio  quickly  reach  an  equilibrium  value  in  the 
Pachappa  sand  but  slowly  decrease  with  time  in  the  Chino  Clay.  The  asymptotic, 
diurnally-averaged  value  of  the  ratio  of  the  local  to  the  bulk  gradient  appears 
to  be  about  3.5  for  Chino  Clay  but  closer  to  2.5  in  the  case  of  Pachappa  sand. 
The  ratio  of  gradients  and  curvature  decrease  during  non-evaporating  periods 
because  large  near-surface  gradients  induced  by  the  evaporation  are  spreading 
downward  due  to  transport  of  soil  water,  and  strong  surface  gradients  are  no 
longer  generated  by  surface  evaporation.  In  fact,  during  non-evaporating 
periods  in  the  Pachappa  sand,  the  bulk  gradient  becomes  comparable  to  that  of 
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the  5  cm  gradient  corresponding  to  a  condition  where  the  vertical  moisture 
gradient  varies  only  slowly  with  depth  and  truncation  errors  are  minimal.  When 


the  potential  evaporation  is  increased  from  .1  mm/hr  to  .3  mm/hr ,  an  appropriate 


geophysical  upperbound,  the  ratio  of  gradients  and  the  curvature  are  higher  but 
behave  in  a  similar  fashion. 

The  smaller  ratio  of  gradients  and  smaller  curvature  in  the  case  of 
Pachappa  sand  reflect  the  greater  penetration  of  the  influence  of  drying  due  to 
larger  diffusivity.  In  the  case  of  constant  diffusivity,  the  penetration  depth 
scale  for  the  influence  of  sinusoidally  varying  evaporation  is  /2D/w  where  u  is 
the  frequency  of  variation.  For  a  diurnal  time  scale  and  variations  of 
hydraulic  diffusivity  for  Chino  Clay,  ranging  from  10“^cm2s~^  to 
10"3cm2s_l,  the  penetration  depth  scale  ranges  from  1.7  cm  to  5.2  cm.  For 
variation  of  hydraulic  diffusivity  in  Pachappa  sand,  ranging  from  3xl0“*s~l, 
to  10-1cm2s“l,  the  diurnal  penetration  depth  scale  ranges  from  2.8  cm  to 
52  cm. 

These  rough  calculations  and  the  behavior  of  the  ratio  of  gradients  (Fig. 

4)  suggest  that  truncation  errors  in  the  computed  flux  for  clay  soils  are  large 
unless  the  upper  layer  is  made  quite  thin.  This  expectation  is  demonstrated  in 
Section  5. 

4.  Direct  evaporation 

In  meteorological  models,  the  surface  moisture  flux  is  usually  parameter¬ 
ized  in  terms  of  a  single  soil  water  content  as  suggested  by  Thornwaite  and 
Mather  (1955)  and  Budyko  (1956)  (see  Eagleson  (1982)  for  a  recent  survey  of  such 
models).  As  the  volumetric  water  content  decreases  below  the  saturation  value, 
evaporation  is  usually  modelled  to  continue  at  the  potential  rate  until  the 
water  content  decrease  below  a  critical  or  reference  value.  Then  the  modelled 
evaporation  is  assumed  to  decrease  linearly  with  decreasing  water  content 
vanishing  at  some  small  nonzero  value  of  the  water  content.  Such  models  appear 
to  be  in  qualitative  agreement  with  observations  collected  in  Fig.  5. 
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Examples  of  observed  evaporation  scaled  by  the  potential  value  as  a 
function  of  soil  water  content.  Thick  lines  represent  cases  where 
transpiration  was  thought  to  be  of  minimal  importance.  Individual 
curves  are:  (M) ,  5  cm  layer  of  50-60%  sand  covered  by  little  vege¬ 
tation  (Marsh  et  al . ,  1981);  (DA),  5  cm  layer  of  sandy  loam  covered 
by  ryegrass  (Davies  and  Allen,  1973);  (R)  1.5  m  layer  of  sandy 
loam,  vegetation  recently  burned  (Rouse  et  al.,  1977);  (D) ,  1  or  2 
cm  layer  of  soil,  vegetation  recently  burned  (Barton,  1979);  (W) ,  5 
cm  layer  of  soil  covered  by  dry  rangeland  grass  (Williams  et  al., 
1978);  (MN) ,  10-20  cm  layer  of  heavy  clay  covered  by  shallow  rooted 
grass  (Mukammal  and  Neuman,  1977);  (BD) ,  25  cm  layer  of  sandy  loam 
with  soybean  crop  (Bailey  and  Davies,  1981);  (SC),  37  cm  layer  of 
sandy  clay  loam  with  wheat  crop  (Seaton  et  al.,  1977);  (S) ,  37  cm 
layer  of  sandy  loam  with  wheat  crop  (Seaton  et  al.,  1977). 
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The  dependence  o£  evaporation  on  water  content  as  determined  from  B  is 
illustrated  in  Fig.  6.  The  relationship  between  evaporation  and  averaged  soil 
water  content  is  very  sensitive  to  the  averaging  depth.  The  critical  water 
content  where  evaporation  begins  to  decrease  with  decreasing  soil  water  occurs 
at  a  higher  moisture  value  when  a  deeper  layer  is  considered  (Fig.  6).  That  isr 
the  drying  of  a  surface  sub-layer  decreases  the  evaporation  even  though  the 
average  water  content  of  the  entire  layer  may  suffer  only  a  slight  decrease. 
These  results  imply  that  attempts  to  develop  models  of  evaporation  based  on 
observations  of  volumetric  soil  water  must  take  into  account  any  differences 
between  model  levels  and  soil  observational  levels. 

The  relationship  of  evaporation  to  layer-averaged  soil  water  content  is 
more  complicated  in  the  case  of  on-off  evaporation.  On  succeeding  days  of 
evaporation,  the  evaporation  first  starts  to  decrease  at  a  higher  soil  water 
content  compared  to  the  previous  day.  This  hysteresis  is  due  to  decreasing 
moisture  at  lower  levels  and  generally  smaller  and  more  diffuse  gradients. 

For  similar  reasons  the  relationship  between  surface  evaporation  and 
layer-averaged  soil  moisture  also  depends  on  the  potential  evaporation  rate  *P 
(Fig.  7).  In  particular,  E/Ep,  for  a  given  layer -aver aged  soil  water  content, 
decreases  with  increasing  potential  evaporation.  High  potential  evaporation 
leads  to  strong  gradients  near  the  surface  so  that  the  water  content  at  the 
surface  is  less  for  a  given  layer -aver aged  value.  With  high  potential 
evaporation  rate,  the  soil  has  greater  difficulty  transporting  water 
sufficiently  fast  to  meet  the  atmospheric  demand. 

These  results  demonstrate  some  of  the  disadvantages  of  choosing  the  upper 
soil  layer  to  be  too  thick  when  relating  evaporation  to  soil  moisture;  namely 
that  the  evaporation  becomes  over -sensitive  to  the  layer-averaged  moisture  value 


Li  jlI*  l"  , 


.80 


.60 


.40 


.20 


.  1  » 
I  I 
*  » 
1  I 
1  * 
1  < 
1  I 

I  '  » 

I  I  I 
I  I 
1  * 

I  »  * 

9  i  > 

2  ll 

<  !  i 

s }  I  i 

Q-  ll 

$  II 

O  /I 

*  // 

/  / 


I  5  10  /I00. 


I  5  10 


100 


.40 


Evaporation  as  predicted  by  model  B  for  constant  potential  evapora 
tion  v-  .1  cm/hr  as  a  function  of  the  volumetric  water  content 
averaged  over  layers  of  various  depth  labelled  in  cm. 
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Fig.  7.  Evaporation  as  predicted  by  model  B  for  constant  potential  evapora' 
tion  rates  of  Ep  ■  .1  cm/hr  and  .3  cm/hr,  as  a  function  of  the 
average  volumetric  water  content  9  of  the  upper  5  cm. 
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would  be  required  near  the  surface  where  vertical  gradients  are  large  during 
evaporation.  In  practice,  the  surface  water  content  and  vertical  gradient  roust 
be  defined  over  a  finite  depth  which  is  unfortunately  large  enough  to  cause 
ambiguity  between  evaporation  and  water  content.  With  a  thinner  surface  layer, 
the  measurement  of  water  content  becomes  undefinable  due  to  both  the  large 
vertical  and  horizontal  gradients  which  typically  occur  near  the  surface. 

An  alternate  approach  is  to  employ  an  expression  which  assumes  equilibrium 
between  the  soil  water  content  at  the  surface  and  the  relative  humidity  of  the 
adjacent  air  as  applied  in  McCumber  and  Pielke  (1981)  and  Camillo  et  al. 

(1983).  The  predicted  relative  humidity  allows  determination  of  the  surface 
evaporation  using  a  bulk  aerodynamic  formulation. 

There  is  no  guarantee  that  the  soil  can  meet  the  demand  imposed  by  the  bulk 
aerodynamic  relationship  even  though  this  demand  for  unsaturated  soil  is  less 
than  the  potential  evaporation  predicted  from  the  bulk  aerodynamic  relationship 
with  100%  surface  relative  humidity.  McCumber  and  Pielke  (1981)  proceed  with  an 
iterative  procedure  which  allows  matching  of  near-surface  soil  water  flux  and 
the  evaporation  into  the  atmosphere. 

Another  approach  can  be  constructed  by  assuming  that  evaporation  proceeds 
at  the  potential  rate  until  the  surface  soil  water  content  decreases  to  an 
"air-dry"  value  (Hanks  et  al . ,  1969;  Nimah  and  Hanks,  1973;  Feddes  et  al. , 

1974).  At  the  air-dry  value,  the  water  film  is  so  thin,  perhaps  only  one 
molecule  thick,  that  electrostatic  forces  prevent  evaporation  from  continuing  at 
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the  potential  rate1.  Then  the  evaporation  is  determined  by  the  flux  from  below 
and  consequently  less  than  the  potential  rate. 

This  concept  was  implemented  through  an  iterative  algorithm  in  model  B 
where  the  surface  water  content  and  soil  water  flux  are  iteratively  adjusted. 

In  the  two-layer  model ,  we  proceed  by  replacing  the  iterative  procedure  with  an 
analogous  two  step  approach.  We  first  tentatively  demand  the  evaporation  to  be 
potential  in  which  case 

D(0>[6sfc  '  V/lU1/2)1  +  K(0)  *  Ep  *  (8) 


The  left-hand  side  is  the  finite  difference  estimate  of  the  water  flux.  We  then 
compute  6sfc  from  (8),  where  9SfC  is  that  surface  water  content  required  to 
induce  the  necessary  gradient  to  produce  a  soil  water  flux  equal  to  the  poten¬ 
tial  evaporation  rate.  If  the  resulting  6sfc  is  greater  than  the  air-dry 
value  Bd,  then  the  evaporation  rate  is  equated  to  the  potential  rate.  If 
0afc  is  less  than  the  air-dry  value  %,  in  violation  of  physical 
expectations,  the  potential  evaporation  cannot  be  met.  Then  the  surface 
evaporation  rate  is  controlled  by  the  soil  water  profile  and  the  surface  water 
content  is  equal  to  the  air-dry  value  in  which  case 

E  »  D(8)ted  -  81]/(z1/2)  +  K( 9)  <  Ep  .  (9) 

*With  this  particular  interpretation  of  the  air-dry  soil  water  content,  the 
increase  of  density  of  water  molecules  in  the  first  5-10  molecular  layers  from 
the  soil  particle  is  assumed  to  not  significantly  influence  the  evaporation 
rate.  The  particular  air-dry  value  defined  in  terms  of  a  single  molecular  layer 
can  be  computed  from  the  estimate  of  surface  area  of  the  soil  particles  which  is 
an  increasing  function  of  percent  clay  content. 
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Comparisons  with  B  indicated  that  using  D  and  K  values  averaged  between  the 
surface  and  upper  layer  mid-level  performed  well.  One  could  also  analytically 
integrate  the  flux  relationship  between  the  surface  and  the  mid-level  of  the 
upper  layer.  However,  this  did  not  improve  the  results  of  the  2-layer  model. 

The  air-dry  value  was  chosen  to  be  .05  for  Pachappa  Sand  and  .16  for  Chino  Clay 
as  in  B.  As  in  Feddes  et  al.  (1974),  model  results  were  found  to  be  insensitive 
to  the  choice  of  air-dry  value. 

In  contrast  to  use  of  an  air-dry  value,  the  approach  based  on  an  equilib¬ 
rium  humidity  relationship  leads  to  an  immediate  reduction  of  evaporation  after 
the  surface  soil  water  decreases  below  saturation.  However,  this  reduction  is 
generally  very  slight  until  the  soil  surface  becomes  quite  dry  depending  on  the 
choice  of  coefficient  values  required  for  this  approach.  Then  this  approach 
behaves  similarly  to  the  method  defined  by  (8-9). 

It  must  be  remembered  that  the  surface  soil  water  content  or  surface  matric 
potential  used  in  these  two  evaporation  approaches  is  an  internal  model  para¬ 
meter  with  no  clear  geophysical  analog.  At  the  actual  air-soil  interface,  soil 
water  content  and  relative  humidity  are  characterized  by  extremely  large 
gradients  in  both  the  horizontal  and  vertical  so  that  such  properties  are  diffi¬ 
cult  to  define.  That  is,  the  distribution  of  H2O  molecules  changes  signifi¬ 
cantly  even  on  a  molecular  scale.  However,  the  adjustment  of  the  "surface"  soil 
water  content  or  matric  potential  in  the  above  two  approaches  does  not  signifi¬ 
cantly  violate  the  water  budget  since  these  surface  properties  sure  theoretically 
defined  over  an  infinitesimally  thin  layer. 

The  choice  between  the  above  two  approaches  appears  to  be  less  important 
than  the  choice  of  hydraulic  properties.  In  this  study  we  use  the  latter 
approach  (8-9)  since  it  facilitates  evaluation  of  truncation  errors  by  compari¬ 
son  with  B  and  because  it  is  simpler.  In  fact,  we  will  march  in  time  by  using 
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hydraulic  properties  from  the  previous  time  step  in  which  case  iteration  can  be 
avoided  completely  with  use  of  (8-9).  This  leads  to  substantial  savings  of 
computer  time  but  has  the  disadvantage  of  potentially  large  truncation  errors. 
These  errors  result  from  the  fact  that  the  hydraulic  diffusivity  may  vary  sub¬ 
stantially  over  a  single  time  step.  Comparisons  in  the  next  section  indicate 
that  such  truncation  errors  are  not  serious  with  modest  time  steps. 

5.  Model  comparison 

Clearly,  a  simple  two-layer  model  of  soil  hydrology  must  necessarily  suffer 
significant  truncation  errors  in  the  flux  between  layers  and  could  be  considered 
inadequate  for  many  purposes.  However,  errors  in  the  instantaneous  between- 
layer  flux  generally  translate  into  smaller  errors  in  the  surface  evaporation 
particularly  in  the  case  of  diurnal  variations.  The  evaporation  depends  on  the 
soil  moisture  near  the  surface  which  is  more  sensitive  to  long  term  errors  in 
the  between-layer  flux  and  less  sensitive  to  diurnal  variations  of  the  error  in 
this  flux.  Additional  truncation  errors  occur  in  the  implementation  of  the 
surface  moisture  flux  conditions  (8-9).  While  truncation  errors  in  the  two- 
layer  model  should  dominate  the  difference  in  the  predictions  between  the  two- 
layer  model  and  B,  other  differences  include  the  fact  that  model  B  is  158  cm 
deep  whereas  the  two-layer  model  is  one  meter  deep. 

The  between-layer  flux  and  evaporation  for  the  cases  of  Chino  Clay  and 
Pachappa  sand  are  shown  for  potential  evaporation  specified  to  be  constant 
(Pigs.  8-9)  and  specified  with  on-off  dependence  on  time  (Figs.  10-11).  For 
Pachappa  sand  forced  by  constant  potential  evaporation,  the  two-layer  model 
predicts  the  total  upward  flux  between  layers  to  be  comparable  to  the  5  cm  flux 
in  B  for  the  first  24-hour  period  although  some  phase  error  occurs.  After  the 
first  day,  the  2-layer  flux  decreases  more  slowly  with  time  than  that  in  the 
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Evaporation  rate  (thick  solid)  and  between-layer  flux  (thin  solid) 
for  Pachappa  sand  forced  by  constant  potential  evaporation  of  .1 
cm/hr  for  the  two-layer  model  with  an  upper  layer  of  5  cm.  Also 
shown  are  the  corresponding  evaporation  rate  from  B  (thick  dashed) 
and  the  5  cm  water  flux  (thin  dashed). 
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Evaporation  rate  (thick  solid)  and  between-layer  flux  (thin  solid) 
for  Chino  Clay  forced  by  constant  potential  evaporation  of  .1  cm/hr 
for  the  two-layer  model  with  an  upper  layer  of  1  cm  thickness. 

Also  shown  are  the  corresponding  evaporation  rate  from  B  (thick 
dashed)  and  1  cm  water  flux  (thin  dashed). 


FLUX  (cm/hr) 


r  I  «  J1"  1 1 


71 


o 


< 

tr 

2 

% 

Ul 


Q 

Ul 

_1 

< 

o 

to 


Fig.  10.  Evaporation  rate  (thick  solid)  and  between-layer  flux  (thin  solid) 
for  Pachappa  sand  forced  by  on-off  evaporation  for  the  two-layer 
model  with  a  5  cm  upper  layer.  Also  shown  are  the  corresponding 
evaporation  rate  from  B  (thick  dashed)  and  5  cm  flux  (thin  dashed). 
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Fig.  11.  (a)  Evaporation  rate  (thick  solid)  and  between-layer  flux  (thin  so 
for  Chino  Clay  for  on-off  evaporation  with  a  1  cm  upper  layer. 
Also  shown  are  the  evaporation  rate  from  B  (thick  dashed)  and 
corresponding  1  cm  flux.  After  72  hours  the  two  evaporation  r 
generally  coincide. 
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Fig.  11.  (b)  Same  as  (a)  except  for  a  5  cm  upper  layer  and  5  cm  flux  from  B. 
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high-cesolution  model  B.  This  overestimation  of  upward  water  flux  eventually 
leads  to  an  overestimation  of  the  evaporation  rate  by  about  30%.  The  error  in 
the  total  accumulated  evaporation  is  significantly  less. 

When  the  modelled  Pachappa  sand  is  forced  by  on-off  evaporation,  the  upward 
water  flux  between  the  two  layers  lags  the  5  cm  flux  in  B  by  several  hours 
(Fig.  10).  However,  the  underestimation  of  the  "daytime"  flux  is  approximately 
compensated  by  an  overestimation  of  flux  during  non-evaporation  periods.  As  a 
result,  the  24-hour  between-layer  flux  predicted  by  the  two-layer  model  is  close 
to  that  predicted  by  the  high-resultion  model  B.  Both  models  predict  near 
potential  evaporation  for  the  duration  of  the  numerical  experiment. 

For  the  case  of  modelled  Chino  Clay,  the  depth  of  the  upper  layer  is 
reduced  from  5  cm  to  1  cm  while  the  depth  of  the  lower  layer  is  reduced  from  95 
cm  to  19  cm.  This  change  recognizes  the  small  diffusivity  and  small  diurnal 
penetration  depth  in  clayey  soils  as  discussed  in  Section  3.  Note  that  the 
diffusivity  for  Chino  Clay  is  particularly  small  (Fig.  2).  If  the  original  5  cm 
and  95  cm  thicknesses  are  retained  for  Chino  Clay,  the  two-layer  model 
substantially  underestimates  the  between-layer  flux  and  evaporation  as  is 
illustrated  for  the  case  of  on-off  evaporation  (Fig.  lib).  With  the  reduced 
layer  thickness  (Fig.  11a),  the  underestimation  of  the  between-layer  flux  is 
substantially  less  although  still  significant.  The  evaporation  rate  predicted 
by  the  two-layer  model  and  B  become  almost  identical  after  the  first  two  days. 
The  evaporation  rates  predicted  by  the  two  models  are  also  in  relatively  close 
agreement  for  Chino  Clay  forced  by  constant  potential  evaporation  with  the  same 
1  cm  thickness  for  the  upper  layer  (Fig.  10). 

An  additional  on-off  evaporation  experiment  was  conducted  where  evaporation 
on  the  fourth  day  is  replaced  with  precipitation  at  a  rate  of  .3  cm/hr.  The 


two-layer  model  somewhat  underestimates  the  resulting  downward  flux  of  soil 
water  for  both  soil  types.  Results  for  Pachappa  sand,  which  responds  more 
rapidly  than  Chino  Clay,  are  shown  in  Fig.  12.  Recall  that  in  the  case  of  down¬ 
ward  water  flux,  formulation  (6-7)  automatically  uses  the  hydraulic  diffusivity 
and  conductivity  of  the  upper  layer.  This  substantially  reduces  the  net  trunca¬ 
tion  error. 

The  truncation  errors  revealed  by  the  above  comparisons  can  be  reduced 
further  by  adjusting  the  depth  of  the  upper  layer  although  this  adjustment  is 
dependent  on  situation.  As  the  thickness  of  the  upper  layer  is  decreased,  the 
time  step  must  be  reduced.  However,  for  most  applications,  the  above  truncation 
errors  will  be  smaller  than  errors  resulting  from  the  estimate  of  soil  type  and 
corresponding  dependence  of  diffusivity  on  soil  water  content.  From  the  above 
model  comparisons,  truncation  errors  in  the  between-layer  flux  appear  to  be 
generally  less  than  30%,  particularly  when  averaged  over  a  24-hour  period.  For 
comparison,  inspection  of  Fig.  2  indicates  that  minor  errors  in  the  estimation 
of  soil  type  or  percent  clay  content  would  lead  to  larger  errors  in  the  water 
flux. 

6.  Conclusions 

At  the  onset,  we  concluded  that  at  least  two  soil  layers  are  required  in 
order  to  simultaneously  include  the  dependence  of  evaporation  on  near-surface 
soil  water  content  and  include  the  storage  of  soil  water  which  occurs  in  a 
deeper  layer.  Since  more  than  two  layers  cannot  be  accommodated  in  many  model¬ 
ling  situations,  we  have  evaluated  the  behavior  of  truncation  errors  in  the 
two-layer  model  by  comparing  with  the  high  resolution  model  of  Boersma  et  al. 
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These  comparisons  indicate  that  the  two-layer  model  will  typically  under¬ 
estimate  the  water  flux  within  the  soil  (Section  5)  due  to  truncation  errors. 
However,  the  impact  of  such  truncation  errors  on  interfacial  water  flux  can  be 
minimized  by  choosing  the  upper  layer  to  be  sufficiently  thin  and  by  using  the 
diffusivity  corresponding  to  the  water  content  of  the  layer  from  which  the  water 
flux  originates.  The  reduced  truncation  errors  are  still  significant  and 
increased  resolution  would  be  required  for  many  research  situations.  However, 
the  truncation  errors  are  probably  small  compared  to  errors  in  estimating  the 
behavior  of  the  hydraulic  diffusivity  for  actual  soil  applications  without 
extensive  measurements.  In  particular,  area-averaged  soil  properties  must  be 
estimated  for  use  in  atmospheric  and  hydrological  models.  Since  the  soil 
hydrology  is  very  nonlinear,  it  is  not  clear  how  to  determine  such  an  average 
even  if  sufficient  soil  data  were  available.  Soils  often  exhibit  extreme  hori¬ 
zontal  variations  on  small  scales. 

In  Section  4,  a  surface  moisture  flux  algorithm  was  developed  which  allows 
control  by  the  soil  moisture  profile  near  the  surface  yet  does  not  require 
iteration.  This  procedure  is  found  to  be  preferable  to  the  usual  practice  of 
relating  evaporation  to  the  layer -averaged  soil  moisture  content. 

Since  transpiration  is  related  to  soil  moisture  in  the  root  zone  and  not 
sensitive  to  the  surface  value,  it  could  be  most  simply  included  by  evaluating 
formulations  similar  to  those  surveyed  in  Eagleson  (1982)  (see  also  Fig.  5) 
using  the  layer -aver aged  moisture.  In  such  a  manner,  a  model  of  evapotranspira- 
tion  could  be  constructed  from  the  present  soil  model  which  would  be  suitable 
for  many  modelling  situations  where  simplicity  is  required. 
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7.  Application  to  the  AFGL  model 

It  is  our  opinion  that  the  physics  o£  the  one-layer  model  is  sufficiently 
inadequate  to  be  of  use  in  the  AFGL  general  circulation  model.  The  two-layer 
model  is  able  to  include  the  most  important  physics,  namely  that  the  evaporation 
is  related  to  the  near-surface  soil  moisture  while  storage  of  water  is  related 
to  the  soil  moisture  in  a  deeper  layer. 

Unfortunately,  the  two-layer  model  suffers  significant  truncation  errors. 
However,  these  truncation  errors  are  not  likely  to  be  larger  than  uncertainties 
in  the  variation  of  soil  type  both  geographically  and  with  respect  to  depth. 

The  initial  soil  model  will  contain  "universal"  properties  which  will  be 
comparable  to  those  of  loam  of  intermediate  clay  content. 

The  role  of  vegetation  and  organic  debris  (dead  plant  material  on  the  soil 
surface)  also  leads  to  significant  uncertainties  in  the  modelling  of  the 
hydrological  budget.  As  a  result  of  this  uncertainty  and  difficulty  in 
estimating  soil  type,  we  have  concluded  that  reducing  truncation  errors  by 
installing  more  than  two  layers  does  not  justify  the  increased  computer  time. 

For  this  reason,  future  research  on  the  soil  models  of  general  circulation 
models  would  best  concentrate  on  parameterizing  vertical  structure  through 
profile  functions  and  should  also  concentrate  on  the  behavior  of  moisture 
transfer  content  in  the  immediate  vicinity  of  the  air-soil  interface.  Inclusion 
of  the  thermodynamics  of  the  two-layer  soil  model  is  also  worth  investigating 
further  since  it  is  not  particularly  complicated  and  could  be  important  in 
special  circumstances. 

For  initial  use  in  the  AFGL  soil  model,  we  have  decided  upon  the  Clapp  and 
Hornberger  (1978)  model  for  hydraulic  diffusivity  and  conductivity.  Which  is  of 


the  form 


Where  Kg  is  the  saturation  hydraulic  conductivity,  \(is  the  saturation  matric 
potential,  9a  the  saturation  volumetric  water  content,  and  b,  a  coefficient 
which  depends  on  soil  type.  If  we  assume  a  "universal  soil  type"  to  have 
average  characteristics  comparable  to  sandy  clay  loam  (28%  clay)  then  b*7 .12, 
Ks=l . 6  x  10-5  cm/min,  0a».42,  and  \J/a*29 .9  cm.  *nie  categorization  of 
soil  type  according  to  grid  point  location  is  a  possible  future  refinement. 

The  addition  of  root  uptake  and  transpiration  is  carried  out  in  Chapter  IV, 


Sections  2-3. 
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IV.  Canopy-Transpiration  Budget 


Abstract 

This  chapter  summarizes  extensive  literature  surveys  of  models  of  the 
canopy  water  budget  and  formulations  of  the  relationship  between  actual  evapo- 
transpiration  and  potential  evaporation.  A  new  model  for  the  canopy  water 
budget  is  developed  which  is  simple,  consistent  with  other  parts  of  the  model 
and  avoids  certain  asymptotic  difficulties  of  previous  models. 

The  dependence  of  transpiration  on  potential  evaporation  and  soil  moisture 
deficit  is  formulated  with  a  plant  coefficient  and  a  functional  dependence  on 
the  soil  volumetric  water  content.  This  formulation  is  based  on  a  large  number 
of  observational  studies  reported  here  in  Chapter  IV  and  also  in  Chapter  III. 
Interdependence  between  transpiration,  direct  evaporation  from  the  soil,  and 
evaporation  of  intercepted  water  from  the  canopy  is  formulated  in  a  manner  to 
insure  self-consistency  and  certain  asymptotic  conditions  for  the  total  evapo- 
transpiration. 


1.  Canopy  water  budget 

a.  Introduction 

The  amount  of  water  that  is  captured  and  retained  by  vegetation  and  evapo¬ 
rated  back  to  the  atmosphere  without  adding  moisture  to  the  mineral  soil  is 
defined  as  the  inteception  loss.  This  definition  excludes  water  that  is 
temporarily  captured  by  the  vegetation  before  dripping  to  the  ground  or  flowing 
down  the  stems  and  trunks  to  the  ground.  Transpiration  is  not  included  as 
interception. 

Tree  canopies,  grasslands,  crops  and  decaying  organic  material  on  the 
ground  (litter)  all  intercept  water.  Fig.  1  shows  relative  magnitudes  of  the 
hydrological  cycle  of  a  vegetation- 3oil  system,  including  interception  and 
interception  loss  (shaded  areas).  These  losses  may  be  a  significant  addition  of 
water  vapor  to  the  atmosphere.  Reevaporation  of  intercepted  water  from  a 
conifer  forest  canopy  can  amount  to  20-40%  of  the  total  annual  precipitation 
falling  on  the  stand  and  10-20%  for  a  hardwood  forest  (Zinke,  1967).  In  some 
instances  the  annual  loss  can  range  up  to  70%  (Pearce  et  al^. ,  1980).  For 
individual  light  precipitation  events,  a  dense  canopy  could  intercept  nearly 
100%  of  the  rain  or  snow.  Litter  interception  losses,  not  shown,  would  normally 
be  a  small  quantity,  about  half  the  soil  evaporation  or  less. 

When  water  is  retained  on  a  leaf,  transpiration  is  suppressed.  It  has  been 
suggested  that  evaporation  of  intercepted  water  exactly  compensates  the 
suppressed  transpiration  such  that  there  is  no  net  loss  of  water  (Leyton  and 
Carlisle,  1959;  Penman,  1963).  If  this  is  the  case,  then  water  could  be  evapo¬ 
rated  from  the  wetted  portion  of  the  canopy  and  transpired  from  the  dry  portion 
of  the  canopy  at  the  same  rate,  regardless  of  the  total  wetness  of  the  canopy. 
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Pig.  1.  The  hydrological  cycle  in  a  vegetation-soil  system;  components 

approximately  to  scale  for  annual  rainfall  100  mm  and  evaporation 
500  mm. 


However,  recent  research  has  indicated  that  the  rate  of  evaporation  of  inter¬ 
cepted  water  from  a  canopy  could  be  three  times  or  more  greater  than  the  trans¬ 
piration  rate  (e.g.,  Thorud,  1967;  Waggoner  et  al . ,  1969;  Stewart,  1977;  Singh 
and  Szeicz,  1979),  thus  the  total  evaporation  from  a  canopy  could  exceed  the 
potential  evaporation.  Potential  evaporation  is  a  reference  evaporation  which 
is  modified  by  plant  type  and  maturity,  time  of  year,  and  soil  moisture  content 
to  determine  the  actual  evaporation  of  the  plant.  For  short  vegetation,  such  as 
grasses,  transpiration  rates  of  dry  surfaces  are  large  enough  to  equal 
evaporation  of  intercepted  water,  such  that  there  would  be  nearly  no  net 
addition  of  water  to  the  atmosphere. 


In  this  section  we  will  attempt  to  determine  the  canopy  capacity  of  various 
types  of  vegetation  and  then  proceed  with  a  critical  survey  of  models  available 
with  the  purpose  of  defining  a  model  to  determine  the  amount  of  water  evaporated 
from  a  canopy  back  to  the  atmosphere. 

b.  Canopy  capacity 

The  amount  of  water  retained  on  a  canopy  after  rainfall  and  drip  drainage 

have  ceased  in  the  absence  of  wind  is  defined  as  the  canopy  capacity,  S.  This 

value,  expressed  either  as  depth  of  water  or  depth  of  water  per  projected  unit 

area  of  the  canopy,  is  dependent  on  meteorological  and  climatic  conditions, 

species,  and  leaf  or  needle  characteristics.  Zinke  (1967)  provides  an  extensive 

review  of  capacity  values  for  various  American  species.  Table  1  presents 

Table  1.  Average  canopy  capacity  values,  S,  for  various  vegetation  types,  n  is 
the  number  of  different  values  used  to  determine  S. 


Type  (n) 


S  (mm) 


Reference (s) 


Conifers 

Pine  (8) 

0.95 

Zinke,  1967;  Gash  and  Morton,  1978 

Spruce  (2) 

2.00 

Leyton  et  al. ,  1967 ;  Hancock  and 
Crowther,  1979 

Douglas  Fir 

(1) 

1.20 

Robins,  1974 

Mixture  (1) 

2.00 

Bringfelt  and  Harsmer,  1974 

Hardwoods 

Eastern  No. 

Am. 

Hardwoods  0.30-1.60 

Helvey  and  Patric,  1965 

Mixture 

Leafy  (4) 

0.90 

Zinke,  1967;  Leyton  et  al.,  1967; 

Leafless 

(4) 

0.50 

Thompson,  1972. 

Shrubs  (5) 

1.20 

Zinke,  1967 

Grasses  (2) 

1.30 

Zinke,  1967 

Tropical  Highland 

Forest  (1)  1.00 

Jackson,  1975 

average  canopy  capacity  values  for  rainfall  for  various  vegetation  types.  These 
values  should  be  used  with  caution  since  they  are  averages  and  deviations  should 
be  expected  according  to  genus,  climatic  conditions,  rainfall  duration  and 
intensity,  and  dislodging  forces  such  as  wind. 

Zinke  (1967)  speculates  that  capacity  values  of  1.3  ran  for  rain  for  most 
grasses,  shrubs  and  trees  would  be  a  reasonable  average.  Snowfall  required  to 
attain  canopy  capacity  ranged  from  2.3  ran  to  9.1  ran  (for  4  cases)  for  conifers, 
with  3.8  ran  reported  as  a  reasonable  average.  Snowfall  on  hardwoods  is  not 
reported,  but  should  be  minimal. 

c.  Interception  models 
Regression  models 

Horton  (1919)  introduced  a  model  for  total-storm  interception  of  rainfall 
which  fills  the  canopy  to  capacity  according  to 

I  =  S  +  kEt  (1) 

r 

where  I  is  the  interception  in  water  depth  on  the  projected  crown  area,  S  the 
canopy  capacity  in  the  same  units,  E  the  evaporation  rate  in  water  depth  per 
unit  time,  and  tr  the  duration  of  the  storm.  The  ratio  of  vegetal  surface  to 
its  ground  projection,  k,  is  analogous  to  leaf  area  index  as  defined  by  Ross 
(1975)  and  generally  exceeds  1.0.  This  equation  states  that  the  total  inter¬ 
ception  during  a  storm  is  equal  to  the  canopy  capacity  plus  evaporative  losses; 
i.e.,  evaporation  increases  the  amount  intercepted  when  the  canopy  is  filled  to 
capacity . 
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Eq.  (1)  assumes  that  the  canopy  fills  to  capacity  during  the  storm.  For  a 
storm  where  the  capacity  is  not  reached,  the  amount  intercepted  does  not  exceed 
the  total  rainfall  and  could  be  expressed  as 

I  -  Ptr  (2) 

where  P  is  tlj  precipitation  rate. 

Horton  also  suggests  relating  total  interception  to  total  gross 
precipitation  P^,  according  to 

I  -  a  +  b  Pg  (3) 


where  a  and  b  are  constants  for  individual  species.  Hydrologists  have  used  this 
method  quite  often  to  determine  the  hydrology  of  a  catchment  or  basin, 
determining  a  and  b  through  regression  analysis .  The  disadvantages  of  this 
approach  are  discussed  later. 

Some  of  the  models  presented  in  the  literature  are  shown  in  Table  2. 


Table  2.  Regression  models  relating  interception,  I,  to  precipitation  charac¬ 
teristics.  Pg  is  the  gross  rainfall,  D  duration  and  R  intensity. 


Vegetation 

Model 

Reference 

Northern  Hardwoods  East  of 

100th  meridian  (average 
over  many  models) 

Growing  Season 

1-0.083  Pg  +  0.036 

Helvey  and  Patric, 

Dormant  Season 

1-0.059  Pg  +  0.020 

1965 

Old-growth  Douglas  Fir, 

Oregon 

1-0.17  Pg  +  0.046 

Rothacher,  1963 

Eastern  White  Pine 

(western  No.  Carolina) 

1-0.14  Pg  +  0.06 

Helvey,  1967 

Ponderosa  Pine  Laws. 

1-0.06  Pg  +  0.11 

Rowe  &  Hendrix, 

for  Pg  >  0.5  in. 

1951 

Tropical  Forest 

1-0.85  +  0.542  In  Pg 

Jackson,  1975 

(east  Africa) 

1-0.764  +  0.521  In  D 

+  0.50  In  R  (Pg  in  n*n) 


Zink<j's  (1967)  review  paper  of  U.S.  vegetation  presents  additional  models  and 
interception  amounts  as  percentage  of  gross  precipitation.  Molchanov  (1960) 
presents  results  for  eastern  Europe. 

Jackson  (1975)  related  interception  to  various  precipitation  characteris¬ 
tics  using  linear,  quadratic  and  logarithmic  models.  He  found  that  as  little  as 
20%  and  as  much  as  60%  of  the  variance  of  the  interception  observations  could  be 
explained  by  the  various  models.  The  models  in  Table  2  are  two  that  gave  satis¬ 
factory  results,  explaining  46%  (relating  1  to  Pg)  and  56%  (relating  1  to  D 
and  R)  of  the  variance. 

The  regression  models,  while  simple  to  use,  suffer  some  serious  disadvan¬ 
tages.  First,  all  do  not  apply  until  the  precipitation  exceeds  a  threshold 
value.  In  the  case  of  linear  models,  Pg  must  exceed  a/(l-b)  where  a  is  the 
intercept  and  b  the  slope.  Secondly,  the  regression  coefficients  are  usually 
determined  for  a  specific  stand  of  trees  and  probably  do  not  apply  to  the  same 
species  in  another  climate  or  stand  with  different  characteristics.  Finally, 
the  models  provide  little  insight  into  the  physics  of  the  interception  process. 

Aston  model 

Aston  (1979)  proposed  a  model  for  a  single-storm  event  which  fills  the 
canopy  to  capacity.  The  gross  interception  is  expressed  as 

I  *  s{l  -  exp l (1  -  p)  Pg/Sj}  +  kEtc  (4) 

where  p  is  the  throughfall  fraction  or  fraction  of  the  canopy  which  allows 
precipitation  to  pass  through  to  the  ground,  i.e.,  gaps  in  the  canopy,  and  the 
other  variables  are  as  defined  previously.  The  last  term  on  the  right  is  a 
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general  evaporation  term,  identical  to  the  term  used  by  Horton  (1919)  and 
discussed  earlier.  No  attempt  by  Aston  is  made  to  further  define  this  term. 

Gash  model 

Gash  (1979)  developed  a  model  that  computed  the  total  interception  over  a 
number  of  storms  with  sufficient  time  between  storms  to  allow  the  canopy  and 
trunk  to  dry  completely.  In  addition  he  assumed  meteorological  conditions  from 
storm-to-storm  are  similar  so  that  mean  precipitation  and  evaporation  apply  to 
all  storms,  that  there  is  no  drainage  from  the  canopy  when  C  <  S  and  for  C  >  S, 

C  ■  S  within  20-30  minutes  after  rainfall  ceases,  independent  of  C  at  the  time 
rainfall  ceases  where  C  is  the  actual  canopy  water  content.  The  model  for  total 
interception  loss  over  n-Hn  storms  is 

I  II 

n-Hn  n  m 

I  I.  -  <l-p-p*>  l  P'.  +  (1-P-Pt>  S  Pa 

3-1  3  *  3-1  3  3-1  3  (5) 

III  IV  V 

n  n-Hn-q 

+<E/R)  I  (P  -P*  )  +  qS  +  p..  E  P^ 

j-1  3  3  j-1  3 

where  P  is  the  gross  rainfall,  P'  the  rainfall  necessary  to  fill  the  canopy  to 
capacity,  p  the  free  throughfall  coefficient  for  the  forest  stand,  Pt  the 
water  diverted  and  retained  on  stems  and  trunk,  E  and  R  the  average  evaporation 
and  rainfall,  respectively,  computed  from  those  hours  when  the  canopy  is 
saturated.  is  the  stem  and  trunk  capacity,  n  the  number  of  storms  suffi¬ 
cient  to  fill  the  canopy  capacity,  m  the  number  of  storms  not  sufficient  to  fill 
the  canopy  and  q  the  number  of  storms  of  the  total,  n-Hn,  that  fill  the  trunk  to 
capacity.  Terms  I  and  II  are  contributions  to  total  interception  from  storms 
that  fill  and  do  not  fill  the  canopy  capacity,  respectively,  term  III  in  the 


additional  amount  intercepted  due  to  evaporation  when  the  canopy  is  full  and 
terms  IV  and  V  are  contributions  from  trunk  and  stems,  respectively.  The 
evaporation  from  each  storm  is  estimated  from  the  Monteith  (1965)  equation. 

The  amount  of  water  required  to  attain  canopy  capacity  is  given  by 

P*  =  S(-R/E)  In  { 1- (E/R)  (l-p-p,.)”1}*  (6) 

This  model  is  the  most  complete  of  the  interception  models,  separating  and 
computing  interception  losses  resulting  from  evaporation  during  rainfall  and 
from  canopy  characteristics.  A  drawback  is  the  need  to  know  values  of  p,  pt 
and  St,  for  which  little  or  no  observations  exist  for  various  species.  Use  of 
this  model  would  then  require  a  guess  at  these  values.  Also  long  term 
observations  to  compute  E  and  possible  R  for  forest  stands  may  not  exist  and 
empirical  values  of  E  and  R  may  be  necessary. 

Water  balance  models 

The  amount  of  water  on  a  canopy  is  charged  by  gross  precipitation  and 
depleted  by  evaporation  and  drip  drainage  (or  blowoff  in  the  case  of  snow).  In 
the  previous  section  gross  amounts  of  intercepted  water  per  storm  period  were 
determined,  primarily  by  relating  interception  to  gross  precipitation.  This 
section  will  present  models  that  determine  changes  in  time  of  water  storage  on  a 
canopy  by  considering  drip  and  evaporation  from  the  vegetal  surfaces.  Most  of 
the  models  consider  only  evaporation  of  intercepted  water  and  not  transpiration 
since  they  ate  concerned  with  the  hydrology  of  a  watershed  rather  than  the 
amount  of  water  evaporated  back  to  the  atmosphere.  These  models  could  be 
extended  to  include  transpiration. 


Rutter  model 

Rutter  et^  al.  (1971)  proposed  a  single-layer  canopy  model  to  describe 
interception  loss  in  pine  canopies  by  calculating  a  running  balance  in  time  of 
rainfall,  evaporation,  throughfall  which  includes  drip  drainage,  and  changes  in 
canopy  storage. 

The  change  in  canopy  storage  of  water,  AC,  for  any  period  within  a  storm, 
is  given  by 


AC  ■  ER  -  ET  -  EE 


where  ER  is  total  rainfall,  ET  total  throughfall  and  EE  total  evaporation. 
These  values  are  expressed  as  equivalent  depth  of  water.  Total  throughfall  is 
the  sum  of  total  drip  drainage,  ED,  and  "free"  throughfall  or  rainfall  through 
gaps  in  the  canopy 


ET  *  ED  +  pER 


where  p  is  proportion  of  rain  falling  through  the  canopy.  Substituting  (8)  into 
(7)  yields 


AC  -  (1-p)  ER  -  ED  -  EE. 


Drip  is  assumed  not  to  begin  until  the  amount  of  water  on  the  canopy,  C, 
exceeds  the  capacity  of  the  canopy,  S.  After  C  >  S  then  the  drip  rate  D  (in 
mm/ unit  time)  is  given  by 

D  *  K  exp(bC) 

where  b  is  an  empirically  derived  constant  and  K  *  exp (In (0.002)-  bS) .  When  the 
canopy  water  is  reduced  to  0,  this  parameterization  incorrectly  continues  to 
produce  a  finite  drip,  at  a  rate  equal  to  K.  Improvements  to  this  are  discussed 
later  in  other  models. 

The  evaporation  rate  E  (in  nm/unit  time)  is  assumed  to  be 
(e  if  C  >  S 

E  J  p 

)E  •  C  if  C  <  S 

(  P  S 

where  Ep  is  the  potential  or  reference  evaporation  rate.  This  is  a  crude 
assumption  since  the  canopy  is  assumed  to  be  one  layer  and  the  same  amount  of 
water  C  could  be  distributed  differently,  thus  changing  the  evaporation.  A 
possible  modification  to  this  is  presented  in  Section  3  of  this  Chapter. 

In  both  of  the  above  we  note  that  C  could  exceed  S.  Rutter  neither 
explains  how  this  is  possible  or  by  how  much  C  should  be  allowed  to  exceed  S. 
Physically,  some  excess  water  could  temporarily  be  detained  on  the  canopy  before 
dripping  through.  Evidence  of  canopy  water  exceeding  capacity  is  the  drip 
occurring  after  rainfall  ceases  and  in  the  absence  of  wind  or  other  dislodging 
forces . 

The  final  form  of  Rutter's  model  for  computing  changes  in  canopy  water 
storage  during  rainfall  is 


VVJ, 


dC 

dt 


c  >  s 


(10) 


(l-p)R  -  K  exp(bC)  -  E 


(l-p)R  -  K  exp(bC)  -  E  •  C  C  <  S 

P  S 


After  rainfall,  (l-p)R  vanishes  and  the  change  in  canopy  water  is  due  to  drip 
and  evaporation.  When  C  >  S,  it  was  assumed  that  there  was  no  drip  (IT  «  pER) ; 
therefore,  the  exponential  drip  term  above  could  be  dropped  for  C  <_  S.  Rutter 
argues,  however,  that  the  drip  rate  is  extremely  small,  on  the  order  of  10~3  mm 
hr"1,  and  unimportant  when  C  <  S,  so  retaining  it  should  cause  little  error. 

This  model  forms  the  basis  for  many  other  models.  It  follows  changes  of 
canopy  water  with  time  by  considering  simplified  physics  of  the  interception 
process.  Its  major  drawback  is  its  treatment  of  evaporation  of  water  from  the 
canopy.  No  distinction  is  made  between  transpiring  and  evaporating  leaf 
surfaces  or  the  distribution  of  canopy  water.  Rutter  et  al .  (1975)  considered  a 
more  generalized  interception  loss  model,  but  did  not  improve  the  evaporation 
computation. 


Calder  model 

Calder  (1977)  made  Rutter's  drip  term  linear  to  eliminate  the  finite  drip 
from  a  dry  canopy.  His  model,  using  the  same  notation,  is 


-  -b.C  +  (l-p)R  -  E 


where  b^ 


b^  for  C  >  S 
b2  for  C  <_  S 


and  E  is  determined  from  the  Penman-Monteith  equation.  Using  optimization 
techniques  Calder  determined  the  optimum  parameters  for  bj,,  b2,  p  and  S. 


Simulated  stocm  results  were  compared  to  actual  storms  for  a  spruce  forest  and 


it  was  found  that  S  ■  1.95mm. 

MANTA  model 

The  MANTA  model  (Sellers  and  Lockwood,  1981)  is  a  multi-layer  crop  model 
that  has  inputs  to  and  outputs  from  four  levels  within  a  vegetal  covering. 
Although  too  complex  for  our  needs,  important  features  of  the  model  are  that 
evaporation  from  the  wetted  portion  of  the  leaf  and  transpiration  from  the  dry 
portion  are  both  considered  and  that  there  is  a  quadratic  relationship  between 
wetted  area  of  a  leaf  layer  and  C/S  for  each  layer.  This  latter  assumption  is 
discussed  in  developing  a  model  for  our  use. 


Massman  model 

Massman  (1980)  developed  a  model  to  predict  changes  in  canopy  storage  which 
included  the  components  found  in  previous  models  but  computed  by  other  methods. 
His  model  is  of  the  form 


dC 

dt 


a(C/S)  B(C/S) 

I(t)  -  D  (S - =i)  -  E  (S— 5 - 

ov  a  .  ’  oK  8  ‘ 


e  -1 


e-1 


(12) 


where  I(t)  is  the  interception  intensity,  D0  and  a  drainage  constants,  and 
Eq  and  8  evaporation  constants,  the  latter  four  depending  on  foliage 
characteristics  and  meteorological  conditions.  Eq  could  be  estimated  by  a 
Penman-like  equation,  D(t)  directly  proportional  to  I(t),  i.e.. 


or  I(t)  assumed  to  be  constant. 

A  major  problem  with  this  method  is  the  lack  of  definition  or 
quantification  of  a  and  g  for  any  vegetation  type.  Also  D(t)  proportional  to 
D0  or  I(t)  depends  on  foliar  wetness  and  whether  or  not  it  is  raining.  As 
with  the  Calder  (1977)  model,  the  problem  of  finite  drip  from  a  dry  canopy  has 
been  eliminated. 

NCAR  GCM 

The  National  Center  for  Atmospheric  Research  (NCAR)  general  circulation 
model  contains  a  routine  to  compute  total  intercepted  water  by  a  canopy  within 
its  planetary  boundary  layer  subroutine  (Dickinson  et  al. ,  1981).  The 
computation  is  very  complex  using  numerous  variables  and  attempts  to  include  all 
factors  in  the  interception  process.  The  model,  keeping  with  notation 
conventions  already  established,  is 

H  -  (l-p)R  -  (=f-Etr)  (14) 

where  Ef  and  Etc  are  the  total  evaporation  rate  including  transpiration  and 
the  suppressed  transpiration  of  the  wetted  portion  of  the  canopy,  respectively. 
The  factor  l~p  is  a  calculated  quantity,  dependent  on  subsurface  temperature, 
and  is  given  by 

Of  summer  Tg2  >  298R 

1-p  =  Of  *  Of  summer  -  {l-FSEAS (Tg2> }AOf  273K  £  Tg2  _<  298K 

Of  summer  -  Aof  Tg2  <  27 3R 

where  Of  summer  and  Aof  are  predetermined  and  tabulated  by  vegetation  type. 

FSEAS (Tg2>  is  a  seasonal  adjustment  of  Aof,  depending  only  on  subsurface 
temperature. 


The  NCAR  model  also  computes  the  maximum  amount  o£  water  a  canopy  can  hold 
wdmax •  This  is  similar  to  canopy  capacity  S  except  that  it  apparently 
includes  temporary  storage  due  to  drip  time  lag.  is  a  vegetation  and 

seasonally  adjusted  foliage  cover  coefficient  given  by 

*  0.02cm  C-  LSAI 
dmax  f 


where  LSAI  is  a  leaf -stem  area  index  defined  by  LSAI  *  LAI  +  SAI  (»  leaf  area 
index  +  stem  area  index).  SAI  is  a  tabulated  value  according  to  vegetation  type 
since  it  does  not  depend  on  season  but  LAI  is  computed  by 


LAI  *  LAI  .  +  FSEAS  (T  „)  (LAI  -LAI  .  ) 

min  g2  max  mm 


f (season,  vegetation  type) 


and  LAImax  and  LAIm^n  are  specified  according  to  vegetation  type. 

Precipitation  is  added  to  the  canopy  if  C  <  Wdmax .  otherwise  C  -  w<jroax 
and  any  excess  is  added  to  the  soil  as  rain  or  snow,  depending  on  the  choice  of 
a  temperature  for  rain  or  snow  determination. 

Ef  represents  the  total  evaporation  of  moisture  from  vegetation  to 
atmosphere  which  is  reduced  by  the  total  suppressed  transpiration,  &tr*  Both 
of  the  above  require  many  computations  and  are  not  presented  here.  An  important 
point  in  the  computations  is  that  each  one  is  adjusted  for  the  proportion  of 
canopy  wetness  through  the  factor  (C/W^max) and  also  adjusted  for 
seasonal  foliage  cover.  This  is  one  step  further  than  the  MANTA  model 


previously  discussed  in  that  seasonal  foliage  cover  is  considered. 


Stanford  Watershed  Model 


A  hydrological  model  currently  used  by  people  interested  in  the  water 
balance  of  a  watershed  is  the  Stanford  Watershed  Model,  SWM-IV  described,  in 
part,  in  Viessman  et  al.  (1977).  In  this  model  all  incoming  precipitation  is 
intercepted  unless  the  precipitation  rate  exceeds  the  interception  rate, 
predefined  according  to  vegetation  type  and  cover,  or  if  the  canopy  capacity  is 
attained.  Evaporation  of  intercepted  water  is  at  the  potential  rate  until  the 
interception  storage  is  depleted. 

The  model  seems  to  assume  that  all  surfaces  of  the  foliage  are  wetted  or 
dried  at  the  same  rate  and  that  canopy  water  storage  does  not  exceed  the  canopy 
capacity. 


d.  Litter  interception 

Litter,  the  dead  and  decaying  organic  material  on  the  forest  floor,  has  not 
received  much  attention.  The  few  results  available  indicate  that  the  amount  of 
water  intercepted  by  litter  is  a  function  of  stand  age.  It  is  estimated  that 
2-5%  of  the  gross  annual  rainfall  is  evaporated  from  litter  (Zinke,  1967).  For 
our  purposes  interception  loss  from  litter  will  be  assumed  negligible.  A  more 
important  function  of  litter  interception  is  to  reduce  direct  evaporation  from 
the  soil  which  may  be  necessary  to  include  in  the  soil  moisture  model. 


2.  Transpiration 


For  the  limited  complexity  which  can  be  accommodated  here,  transpiration  is 
best  related  to  potential  evaporation  by  means  of  a  single  coefficient  or 
function  without  explicit  representation  of  stomatal  and  internal  plant  resis¬ 
tances.  This  coefficient  is  often  referred  to  as  the  crop  coefficient  or  plant 
coefficient.  This  approach  is  the  usual  one  in  routine  irrigation  management 
and  simple  hydrological  modelling.  As  a  result  there  are  many  studies  attempt¬ 
ing  to  determine  the  value  of  the  plant  coefficient;  some  of  these  studies  are 
summarized  in  Table  3. 

Table  3.  Plant  coefficients  (PC)  for  different  seasons  at  different  locations. 

Maximum  refers  to  the  maximum  for  the  growing  season;  growing  season 
is  the  average  for  the  growing  season.  General  refers  to  any  location 
where  modification  by  empirical  coefficients  of  climate  and  geography 
are  necessary  for  application  of  PC  to  a  specific  location.  Necessary 
for  the  determination  of  PC  is  the  way  in  which  PET  and  ET  are 
obtained.  Comb,  refers  to  a  combination  equation,  TB  is  a  thermally 
based  equation,  pan  is  pan  evaporation,  ETB  is  an  energy  and  thermally 
based  equation,  Lys  is  a  lysimeter  measurement  and  BC  refers  to  an 
already  existing  Blaney-Criddle  crop  coefficient  so  that  no  PET 
measurement  was  made. 


Season 

PC 

Location 
where  valid 

PET 

measurement  ET 

First  method  or  measurement 

author.  Date  equation  method 

Agricultural  plants: 

Crop:  Alfalfa 

1. 

Spring 

.93 

Arizona 

Van  Bavel,  1966  comb. 

lys. 

Summer 

.99 

Fall 

.96 

2. 

Maximum 

1.00 

Arizona 

Van  Bavel,  1967  comb. 

lys. 

3. 

Maximum 

1.35 

General 

Hargeaves,  1974  TB 

lys. 

Growing 

1.00 

4. 

Maximum 

.85 

Semi -arid  & 

Holmes,  1959  - 

BC 

Arid  locations 

5. 

Maximum 

1.10 

General 

USDA* ,  1967 

BC 

Yearly  Avg. 

.80 

Table  3.  (continued) 


Crop:  Snap  Beans 

1. 

2.  Sumner 

3 .  Maximum 
Growing  Season 

4 .  Maximum 
Growing  Season 

Crop:  Dry  Beans 

1.  Maximum 
Growing  Season 

2 .  Maximum 
Growing  Season 


Kimberly,  Idaho 

.48  Wisconsin 

1.10  General 
.85 

1.15  General 

1.10  General 
.85 

1.15  General 
.90 


Wright,  1978  oomb 
Black,  1970  comb 
USDA,  1967  - 

Hargreaves ,  TB 
1974 


USDA,  1967 


Hargreaves ,  TB 
1974 


Table  3.  (continued) 


Crop:  Soy  Beans 


1. 

Maximum 

1.05 

General 

USD*,  1967 

BC 

Growing  Season 

.65 

2. 

Maximum 

.73 

India 

Vankatachari,  - 

BC 

Growing  Season 

1978 

3. 

Maximum 

1.15 

General 

Hargreaves,  TB 

lys. 

Growing  Season 

1974 

4. 

Maximum 

1.10 

Iowa 

Stanley,  1978  pan 

neutron-probe 

Growing  Season 

.70 

gravimetric 

Cco] 

a:  Corn 

1. 

Summer 

1.15 

Ohio 

Parmele,  1974  comb. 

lys. 

Late  Summer 

1.01 

2. 

Maximum 

1.25 

Israel 

Lomas,  1974  comb. 

lys. 

Growing  Season 

.80 

3. 

Maximum 

1.15 

Minnesota 

Dylla,  1980  BTB 

lys. 

Growing  Season 

.70 

4. 

Maximum 

.89 

India 

Vankatachari,  - 

BC 

Growing  Season 

.71 

1978 

5. 

Maximum 

1.15 

General 

Hargreaves,  TB 

lys. 

Growing  Season 

1974 

6. 

Maximum 

1.08 

General 

US DA,  1967 

BC 

Growing  Season 

.85 

Crop:  Cotton 

1. 

Maximum 

.70 

Israel 

Fuchs,  1963 

BC 

2. 

Maximum 

.81 

India 

Vankatachari,  - 

BC 

Growing  Season 

.60 

1978 

3. 

Maximum 

1.02 

General 

USDA,  1967 

BC 

Growing  Season 

.60 

4. 

Maximum 

1.15 

General 

Hargreaves,  TB 

lys. 

Growing  Season 

.90 

1974 

101 


MSI 


s 


g 


vfl 


& 


Table  3.  (continued) 

Crop:  Grain  (wheat) 


1. 

Maximum 

.81 

India 

Vankatachar i. 

— 

BC 

Growing  Season 

.72 

1978 

2. 

Maximum 

.95 

Eastern  Oregon 

Bates,  1982 

pan 

neutron 

Growing  Season 

.70 

probe 

3. 

Maximum 

1.30 

General 

USDA,  1967 

— 

BC 

Growing  Season 

.75 

4. 

Maximum 

1.15 

General 

Hargreaves , 

TB 

lys. 

Growing  Season 

.90 

1974 

Cro| 

3:  Grain  (Barley) 

1. 

Growing  Season 

.76 

Denmark 

Kristensen, 

Penman 

neutron 

1974 

probe 

Crop:  Grain  (Sorghum) 

1. 

Maximum 

.72 

India 

Vankatachar i , 

— 

BC 

Growing  Season 

.56 

1978 

2. 

Maximum 

1.08 

General 

USDA,  1967 

BC 

Growing  Season 

.70 

3. 

Maximum 

1.15 

General 

Hargreaves, 

— 

lys. 

Growing  Season 

.90 

1974 

Crop:  Grass  (short. 

green) 

1. 

Summer 

.80 

SE  England 

Thom,  1977 

Penman 

? 

Winter 

.60 

2. 

Growing  Season 

.73 

Denmark 

Kristensen, 

Penman 

neutron 

1974 

probe 

3. 

Maximum 

1.00 

General 

Hargreaves, 

TB 

lys. 

Growing  Season 

1.00 

1974 

4. 

Spring 

.75 

Wyoming 

O'Neill,  1979 

— 

BC 

Summer 

Fall 


■COiv.  v.-wV- 


Table  3.  (continued) 


i 

Crop:  Grass  (long)  Zj 


1.  Growing  Season 

Crop:  Grass  (Pasture) 

.74 

Denmark 

Kristensen, 

1974 

Penman 

neutron  y 

probe  'j 

V 

V 
u 

*  U 
*] 

1.  Summer 

.90 

Australia 

Shepherd,  1972  comb. 

iys.  ;j 

2 .  Maximum 

1.50 

Netherlands 

Strieker, 

bulk 

energy 

Growing  Season 

.85 

1978 

aerodynamic 

balance 

3.  Maximum 

Crop:  Potatoes 

1.15 

General 

Hargreaves, 

1974 

TB 

lys.  £ 

t 

1.  Maximum 

Growing  Season 

1.38 

.90 

General 

USDA,  1967 

— 

BC  5 

C 

2 .  Summer 

.90 

Australia 

Shepherd,  1972  comb. 

lys. 

3.  Maximum 

Growing  Season 

.90 

.50 

Eastern  Oregon 

Bates,  1982 

pan 

neutron  f4 

probe  ? 

4.  Maximum 

Growing  Season 

Crop:  Rice 

1.15 

.90 

General 

Hargreaves, 

1974 

TB 

lys .  \ 

. 

S' 

1.  Maximum 

1.20 

Arid  6  Semi-Arid 
Conditions 

Holmes,  1959 

— 

BC  j 

2.  Sumner 

Crop:  Sugar  Cane 

1.02 

California/Davis 

Lourence, 

1971 

reference 

crop 

energy  R 

balance  9 

V 

3 

V 

V 

1.  Maximum 

Growing  Season 

1.25 

1.00 

General 

Hargreaves, 

1974 

TB 

lys.  ^ 

2.  Growing  Season 


85  General 


USDA,  1967 


BC 


Table  3.  (continued) 

Croi 

a:  Sugar  Beets 

1. 

Growing  Season 

.82 

Denmark 

Kr istensen. 

Penman 

neutron 

1974 

probe 

2. 

Maximum 

1.15 

General 

Hargreaves, 

TB 

lys. 

• 

Growing  Season 

.90 

1974 

3. 

Maximum 

1.25 

General 

USDA,  1967 

— - 

BC 

Growing  Season 

.90 

Crop:  Field  and  Oil  Crops 

(Flax,  peanuts 

,  safflower,  tomatoes) 

1. 

Maximum 

1.15 

General 

Hargreaves, 

TB 

lys. 

Growing  Season 

.90 

1974 

Crop:  Citrus  Fruits 

(Orange,  lemon  grapefruit) 

1. 

Maximum 

.75 

General 

Hargreaves, 

TB 

lys. 

Growing  Season 

.75 

1974 

2. 

Maximum 

.72 

General 

USDA,  1967 

— 

BC 

Yearly  Avg. 

.70 

Crop:  Deciduous  Fruits  (peaches,  plums. 

walnuts) 

1. 

Maximum 

1.10 

General 

Hargreaves, 

TB 

lys. 

Growing  Season 

.85 

1974 

2. 

Maximum 

1.00 

General 

USDA,  1967 

— 

BC 

Yearly  Avg. 

.50 

Croj 

3:  Grapes 

1. 

Maximum 

.75 

General 

Hargreaves, 

TB 

lys. 

Growing  Season 

.60 

1974 

2. 

Maximum 

.82 

General 

USDA,  1967 

— 

BC 

Yearly  Avg. 

.50 

«-  *  -  A  .V  .V  J 

Table  3.  (continued) 


Crop;  Sumner  Vegetables 


Maximum 

Growing 

Season 

1.15 

.85 

General 

Hargreaves, 

1974 

TB 

lys 

Maximum 

Growing 

Season 

.85 

General 

USDA,  1967 

— 

BC 

Mon-agricultural  plants; 


Plant:  Spruce, 

Pine,  Fir 

Spring 

1.00  France 

Assenac,  1972  ? 

? 

Summer 

1.00 

Fall 

1.00 

Plant:  Douglas  Fir 

Annual  Range  .64  to  France  Assenac,  1972  ?  ? 

1.00 

Plant:  Prairie  grass 


Summer 

.90  Canadian  Great 

Nkemdirim, 

Penman 

lys/ 

Plains 

1973 

neutron 

probe 

Table  3.  (continued)  Seasonal  plant  coefficients  (PC)  for  irrigated 
crops  from  (JSDA  (1967). 


Length  of  Normal  Growing  Plant 


Crop 

Season  of  Period  (1) 

Coefficient  (2) 

Alfalfa 

Between  frosts 

0.80 

to 

0.90 

Bananas 

Full  year 

.80 

to 

1.00 

Beans 

3  months 

.60 

to 

.70 

Cocoa 

Full  year 

.70 

to 

.80 

Coffee 

Full  year 

.70 

to 

.80 

Corn  (Maize) 

4  months 

.75 

to 

.85 

Cotton 

7  months 

.60 

to 

.70 

Dates 

Full  year 

.65 

to 

.80 

Flax 

7  to  8  months 

.70 

to 

.80 

Grains,  small 

3  months 

.75 

to 

.85 

Grain,  sorghums 

4  to  5  months 

.70 

to 

.80 

Oilseeds 

3  to  5  months 

.65 

to 

.75 

Orchard  crops: 

Avocado 

Full  year 

.50 

to 

.55 

Grapefruit 

Full  year 

.55 

to 

.65 

Orange  and  lemon 

Full  year 

.45 

to 

.55 

Walnuts 

Between  frosts 

.60 

to 

.70 

Deciduous 

Between  frosts 

.60 

to 

.70 

Pasture  crops: 

Grass 

Between  frosts 

.75 

to 

.85 

Ladino  whiteclover 

Between  frosts 

.80 

to 

.85 

Potatoes 

3  to  5  months 

.65 

to 

.75 

Rice 

3  to  5  months 

1.00 

to 

1.00 

Soybeans 

140  days 

.65 

to 

.70 

Sugar  beet 

6  months 

.65 

to 

.75 

Sugarcane 

Full  year 

.80 

to 

.90 

Tobacco 

4  months 

.70 

to 

.80 

Tomatoes 

4  months 

.65 

to 

.70 

Truck  crops,  small 

2  to  4  months 

.60 

to 

.70 

Vineyard 

5  to  7  months 

.50 

to 

.60 

(1)  Length  of  season  depends  largely  on  variety  and  time  of  year  when  the  crop 
is  grown.  Annual  crops  grown  during  the  winter  period  may  take  much  longer 
than  if  grown  in  the  summertime. 

(2)  The  lower  values  of  plant  coefficients  for  use  in  the  Blaney-Criddle 
formula  are  for  the  more  humid  areas,  and  the  higher  values  are  for  the 
more  arid  climates. 
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Although  such  studies  are  ubiquitous,  three  major  uncertainties  remain: 

(1)  The  crop  coefficient  is  related  to  a  variety  of  different  definitions  of 
potential  evaporation  which  might  significantly  influence  the  value  of  the 
plant  coefficient. 

(2)  Values  have  been  estimated  primarily  for  commercial  crops.  Information  on 
natural  plant  communities  is  limited. 

(3)  Due  to  the  usual  observational  set  up,  the  plant  coefficient  usually 
includes  both  transpiration  and  direct  evaporation  from  soil.  He  will  want 
to  distinguish  between  transpiration  and  direct  soil  evaporation  because 
the  relationship  of  these  two  processes  to  the  soil  moisture  distribution 
is  different. 

a.  Plant  coefficient 

The  plant  coefficient  is  formally  defined  here  as  the  ratio  of 
transpiration  to  the  potential  evaporation  for  the  case  of  insignificant  soil 
water  deficit.  Unfortunately,  the  plant  coefficient  is  usually  evaluated  from 
measurements  as  the  ratio  of  the  total  evapotranspiration  to  the  potential 
evaporation.  However,  with  significant  vegetative  cover  the  contribution  due  to 
direct  evaporation  will  usually  be  small.  The  plant  coefficient  accounts  for 
the  various  ways  which  plant  density  and  stomatal  behavior  influence  the  actual 
transpiration.  He  will  later  partition  the  influences  of  vegetative  density  and 
stomatal  resistance.  He  leave  the  two  factors  combined  here  since  the 
distinction  between  them  is  usually  not  possible  from  observations. 

Evapotranspiration  may  be  determined  by  direct  or  indirect  methods.  An 
example  of  direct  measurements  of  (w'q')  is  the  computation  from  observations  of 
fluctuations  of  specific  humidity  q  and  vertical  motion  w. 


Evapotranspiration  may  be  estimated  indirectly  by  means  of  the  energy 
budget.  Net  radiation,  soil  heat  flux  and  sensible  heat  flux  are  determined; 
the  residual  term  is  then  equated  to  latent  heat  flux  (e.g.,  Fritschen,  1965). 
The  soil  moisture  budget  method  equates  evapotranspiration  to  precipitation  less 
change  in  soil  moisture  storage,  runoff  at  the  surface  and  infiltration  to  the 
ground  water  zone  beneath  (Strahler,  1975).  Instrumentation  used  in  measuring 
the  change  in  soil  moisture  include  lysimeters  and  other  gravimetric  devices, 
and  the  neutron  probe.  A  lysimeter  is  an  isolated  block  of  soil  in  which  soil 
moisture  can  be  closely  monitored.  In  general,  gravimetric  methods  involve  a 
"before"  and  "after"  soil  moisture  measurement  which  is  then  related  to  evapo-* 
transpiration.  The  neutron  probe  is  an  instrument  that  measures  the  scattering 
of  neutrons  in  the  soil  and  relates  this  to  water  content. 

b.  Measurements 

Fig.  2  is  an  idealized  example  of  the  annual  variation  of  potential  evapo¬ 
transpiration,  actual  evapotranspiration,  direct  evaporation  from  the  soil,  and 
transpiration.  Available  information  on  evaporation  of  intercepted  water  does 
not  allow  construction  of  a  curve  for  seasonal  variation.  Note  that  the  trans¬ 
piration  curve  peaks  more  sharply  in  summer  than  the  soil  evaporation  or  poten¬ 
tial  evapotranspiration  curves.  This  is  due  to  the  large  increase  in  the  plant 
coefficient  mainly  associated  with  increasing  plant  density.  The  sum  of  the 
soil  evaporation  and  transpiration  curves  gives  the  evapotranspiration  curve 
which  also  peaks  sharply.  As  a  result  the  value  of  the  plant  coefficient  peaks 
in  summer  (Fig.  3).  Note  the  significant  decrease  of  the  plant  coefficient  as 
the  plant  approaches  maturity. 

The  plant  coefficients  in  Table  3  represent  typical  values  from  numerous 
studies.  Listed  are  plant  coefficients  for  vegetation  types  for  different  times 


at: 


Winter 


Spring  Summer  Autumn 


Typical  examples  of  annual  curves  of  potential  evapotranspiration 
(PET),  evapotranspiration  (ET),  soil  evaporation  (E),  and  plant 
transpiration  (T) .  Curves  generalized  from  Al-Khafaf  et  al. (1978); 
Bates  et  al.  (1982);  Cuenca  et  al. (1981) ;  Fritschen  (1982); 
Kanemasu  et  al.  (1976);  Nelson  and  Hwang  (1976);  Russell  (1980). 


Winter  Spring  Summer  Autumn 


Typical  example  of  an  annual  plant  coefficient  (PC)  curve.  Dashed 
.line  indicates  modification  due  to  harvest.  Curve  generalized  from 
Bates  et  al.  (1982)  and  Wright  and  Jensen  (1978). 


of  the  year  at  various  geographic  locations  for  the  condition  of  non-limiting 
soil  moisture. 

c.  Plant  coefficient  partitioning 

From  a  physical  point  of  view  it  is  useful  to  separate  the  influences  of 
stomatal  resistance  and  vegetative  density  because  the  latter  leads  to  reduction 
of  direct  evaporation  from  the  soil  due  to  decreased  solar  radiation  and  air 
flow  at  the  soil  surface. 

Therefore ,  we  represent  the  plant  coefficient  as 

E>p/Ep  *  kv  Of 

where  Bj>  is  the  transpiration,  Ep  the  potential  evaporation,  ky  a  factor 
due  to  stomatal  resistance  analagous  to  that  used  in  Eagleson  (1982)  and  Of  a 
factor  representing  vegetative  density  (Deardorff,  1978;  Eagleson,  1982)  which 
will  be  referred  to  as  the  shading  factor. 

The  coefficient  a f,  which  varies  between  zero  and  one,  has  the  advantage 
that  it  also  represents  the  percent  of  solar  radiation  reaching  the  soil  surface 
after  adjusting  for  sun  angle.  The  principle  disadvantage  is  that  measurements 
of  transpiration  have  been  usually  related  to  the  "leaf  area  index."  Some 
examples  are  shown  in  Fig.  4.  Fortunately  these  results  are  of  guidance  here 
since  the  relation  between  shading  factor  and  leaf  area  index  has  been  estimated 
in  several  studies  (Fig.  5). 

The  product  kyOf  typically  varies  between  .5  and  unity  (Table  3). 
Emphasizing  the  few  observations  over  grasses  and  forest,  we  choose  a  typical 
plant  coefficient  to  be  .7.  In  the  present  model  the  reduction  from  unity  is 
attributed  entirely  to  the  shading  factor  while  the  resistance  factor  is 


Relative  evapotranspiration  as  a  function  of  leaf  area  index  from 
Hanson  (1976)  (open  circle,  solid  line);  Kristensen  (1974)  (solid 
circles,  solid  line);  Ritchie  and  Burnett  (1971)  (x's,  solid  line); 
Eagleson  (1978)  for  ky^l.O  (open  circles,  broken  line)  and 
kv=.7  (x's,  broken  line). 


assigned  to  be  unity.  Such  values  can  be  altered  to  accommodate  special 
circumstances  or  global  variations. 

d.  Dependence  on  soil  moisture  deficit 

Some  examples  of  the  dependence  of  evapotranspiration  on  the  soil  water 
content  are  presented  in  Fig.  5  of  Chapter  III.  While  transpiration  appeared  to 
decrease  continuously  with  the  soil  moisture  deficit,  a  sharp  decrease  begins  to 
occur  at  moderately  low  values  of  soil  moisture.  The  transpiration  continues  to 
decrease  rapidly  until  it  vanishes  at  some  nonzero  value  of  moisture.  The 
latter  is  referred  to  as  the  wilting  point. 

As  noted  in  Chapter  III,  the  evapotranspiration  is  often  modelled  to  be 
potential  until  the  soil  moisture  decreases  below  a  reference  or  critical  value, 
then  to  decrease  linearly,  vanishing  at  the  wilting  point.  We  will  model  the 
transpiration  part  of  the  evapotranspiration  in  this  matter  as  is  explicitly 
recorded  in  Section  3  of  this  chapter. 

e.  Values  of  soil  parameters 

Examples  of  values  of  the  important  quantities  necessary  as  inputs  to  soil 
moisture  for  modelling  evapotranspiration  sure  briefly  summarized  in  Table  4. 
Runoff  parameters  will  not  be  considered  at  this  time.  Soil  water  is  expressed 
in  terms  of  0,  the  volumetric  soil  water  content  (percent  water  by  volume). 

The  following  symbolism  will  be  used:  0Cap  "  field  capacity,  6<>rit  ” 
value  above  which  water  is  evaporated  at  the  potential  rate,  6sat  *  soil 
saturation,  6pwp  *  permanent  wilting  point,  =■  depth  of  the  soil  layer  in 
the  model.  The  following  relationship  will  generally  be  true 


0pwp  <  9crit  <  8cap  <  Qsat 


but  may  depend  on  a  particular  author's  own  definition  of  these  terms.  Fiel 
capacity  has  uncertain  physical  meaning  but  is  included  here  due  to  its  wide 
spread  usage. 


Table  4.  Soil  moisture  parameters 


I.  Baier  and  Robertson,  1966 
ecap  *  0*30 
9pwp  »  0-15 

model  soil  depth  26  in. 

number  of  layers  6 

available  water  (6cap  -  6^)  4  in. 


Carson,  1982 

/0 cap  range  10-30  cm 
/®crit  range  5-11.25  cm 

Clapp  and  Hornberger,  1978 

SOIL  TYPE  % 

CLAY 

0sat 

sand 

3 

0.395 

sandy  loam 

9 

0.435 

loam 

19 

0.451 

clay  loam 

34 

0.476 

sandy  clay 

43 

0.426 

clay 

63 

0.482 

(This  is  only  a  part  of  a  much 

larger 

set  of  soil  types) 

I7TJ-  •TTTTTTrre 


114 


IV.  Deardorff,  1977 


Dq  »  10  cm 

Di  =  50  cm  The  surface  layer,  Dq,  is  contained  within  the 

^crit  *  0.30  total  depth,  D^. 

®cap  ”  0.40 

V.  Deardorff,  1978 


Do  *■  10  cm 

dpffp  *  0.10 

0cap  - 

0.40 

Di  m  50  cm 

Ocrit  “  0.30 

VI. 

Manabe,  1969 

Dg  “  1  m 

9crit  “  *11 

®cap  “ 

.15 

VII. 

Marshall  and  Holmes, 

1979 

SOIL  TYPE 

«  CLAY 

e 

_E2E 

0cap 

8sat 

Sand 

3 

0.02 

0.06 

0.4 

Loam 

22 

0.05 

0.29 

0.5 

Clay 

47 

0.20 

0.41 

0.6 

VIII.  McCumber  and  Pielke,  1981 
Dq  ■  1  n,  14  levels 


SOIL  TYPE 

PWP 

Sand 

0.068 

Sandy  Loam 

0.114 

Sandy  Clay 

0.219 

Peat 

0.395 

.*•  .V. 


Dickinson  et  al.,  1981 


Dg  ■  10  cm  overlapping  layers,  where  Dg  is  the  surface 

■  100  cm  layer 

0*pWp  *  0.30  (volume  of  water/volume  of  water  at  saturation) 

®*cap  *  0.60 

where  6*  is  the  volume  of  water  in  the  layer  /volume  of  water  in  the 
layer  of  saturation.  With  this  definition  0_<9*<1  and  9*sat  *  1* 

Clothier  et  al.,  1981 
dSat  “  0.34  for  fine  sand 
®cap  “  0  .33  for  sandy  loam 

Lee,  1980 


SOIL  TYPE 

0 

J&E 

0 

cap 

Sand 

0.025 

0.100 

Fine  Sand 

0.033 

0.116 

Sandy  Loam 

0.050 

0.158 

Loam 

0.100 

0.267 

Silty  Loam 

0.116 

0.283 

Light  Clay  Loam 

0.133 

0.300 

Clay  Loam 

0.150 

0.317 

Heavy  Clay  Loam 

0.175 

0.325 

KII.  Rutter,  1975 


SOIL  TYPE 

6 

pwp 

6 

cap 

Sand 

0.02 

0.09 

Sandy  Loan 

0.11 

0.27 

Loam 

0.13 

0.34 

Silt  Loam 

0.14 

0.38 

Clay  loam 

0.16 

0.30 

Clay 

0.22 

0.39 

Peat 

0.25 

0.55 

XIII.  Black,  1979 

SOIL  TYPE 

8 

P"P. 

8 

_£5E 

GROUND  COVER 

Sandy  Loam 

0.08 

0.215 

70  cm 

Douglas  fir 

Sandy  Loam 

0.11 

0.213 

85  cm 

Douglas  fir 
salal 

3.  Canopy  water  and  transpiration  model 

Normally  evapotranspiration  is  considered  collectively  without 
distinguishing  between  transpiration  and  direct  soil  evaporation.  Here  we  want 
to  preserve  the  distinction  since  the  direct  soil  evaporation  is  most 
appropriately  related  to  the  soil  moisture  of  an  upper  thin  layer  while  water 
for  transpiration  originates  more  from  the  deeper  root  zone.  The  total 
evaporation  can  be  written  as 

E  -  E„.  +  E_  +  E  *F»E  (15) 

dir  T  C  p 

where  %ir  is  the  direct  evaporation  from  the  soil ,  is  the  transpiration , 

E£  is  the  evaporation  of  precipitation  intercepted  by  the  canopy,  Ep  is  the 
potential  evaporation  and  F  is  a  function  whose  components  are  discussed  below 
and  summarized  in  Chapter  VI. 

a.  Reduction  of  direct  evaporation 

Vegetation  reduces  the  direct  evaporation  from  the  soil  by  shading  the 
ground  and  reducing  the  wind  speed  near  the  ground.  The  reduction  of  wind  speed 
can  be  posed  in  terms  of  increased  surface  roughness  parameter  and  increased 
displacement  height.  One  can  further  include  subcanopy  flow  dynamics. 

The  roughness  length,  displacement  height  and  existence  of  special 
subcanopy  dynamics  depend  on  the  height  and  density  of  the  vegetation.  Both  of 
these  effects  axe  presumably  dependent  on  the  leaf  area  index  or  related  shading 
factor  discussed  in  Section  2.  The  reduction  of  solar  radiation  reaching  the 
ground  surface  through  the  vegetation  can  be  expressed  as  a  function  of  the 
shading  factor;  here  we  do  not  explicitly  include  complexities  due  to  varying 
sun  angle. 


GRAVITATIONAL 

SETTLING 


Schematic  illustration  of  geometry  of  hydrology  model 


To  minimize  the  number  of  parameters  we  relate  both  the  influences  of 


shading  and  wind  speed  reduction  to  the  shading  factor  o f  according  to  the 
format 


E 


dir 


E  ..  (l-o.) 
soil  f 


(16) 


where  Of  is  the  shading  factor,  Egon  is  the  evaporation  from  the  soil  model 
in  the  absence  of  vegetation  as  modelled  in  Chapter  III.  The  linear  dependence 
on  Of  is  based  partly  on  the  expected  dominance  of  the  importance  of  radia- 
tional  forcing  of  direct  soil  evaporation  under  a  partial  canopy  as  assumed  in 
Feddes  et  al.  (1974). 


b.  Transpiration 

Transpiration  is  related  to  the  density  of  vegetation  cover  and  the 
stomatal,  internal  and  root  resistance  of  the  vegetation  and  the  soil  moisture 
content.  Here  we  detail  the  transpiration  corresponding  to  the  two-layer 
version  of  the  model  although  the  model  is  set  up  to  accommodate  up  to  ten 
layers.  These  influences  jure  most  simply  included  with  the  following  format  for 
transpiration  which  is  based  on  discussions  in  Section  2. 


z  ( z  — z  ) 

ET  ■  Ep  Vf  v  ♦  «<v  t1  - ,c/s,n] 


(17) 


g(  8) 


(  1 

0 

> 

)  9-9wilt 

9 

< 

1  **ref  Ewilt 

l  0 

0 

< 

ref 


Here  we  have  assumed  that  the  root  uptake  rate  is  independent  of  depth  within  a 
given  layer.  The  wilting  point,  S^ilt'  where  root  uptake  ceases,  is  assigned 


to  be  .12.  The  plants  do  not  die  if  the  water  content  decreases  to  the  wilting 
point;  transpiration  begins  again  as  soon  as  6  exceeds  Svu.it •  The  parameter 
0ref  is  the  soil  moisture  content  where  the  soil  moisture  deficit  begins  to 
reduce  root  uptake  and  transpiration.  6ref  is  chosen  to  be  .25  which  is 
significantly  below  the  saturation  values  for  most  soil  types.  These  soil 


parameter  values  are  based  on  results  in  Section  2d. 

C  is  the  canopy  water  content  and  S  is  the  canopy  water  capacity,  both  of 
which  are  discussed  in  Section  1.  The  coefficient  ky  is  the  plant  resistance 
factor  chosen  to  be  unity  and  Of  is  specified  to  be  .7.  As  discussed  in 
Section  2,  the  product  of  kyOf  is  similar  to  the  commonly  used  plant 
coefficient.  The  parameter  n  is  chosen  to  be  1/2  to  be  consistent  with  the 
interception  model  discussed  below. 


c.  Interception 

We  have  constructed  a  model  for  the  canopy  water  budget  by  considering 
previous  models  and  any  inconsistencies  as  discussed  in  Section  1  and 
recognizing  the  level  of  simplicity  required  here.  We  propose  the  following 
interception  model 


ofP  -  Ec 
Of(C/S)nEp 


where  C  is  the  canopy  water  storage,  S  the  storage  capacity  of  the  canopy,  here 
chosen  to  be  2  nin,  P  the  precipitation  rate,  Ep  the  potential  evaporation  rate 
and  Of  the  shading  factor.  This  interception  model  is  similar  to  that  of 
Rutter  et  al.  (1971)  except  that: 


1)  the  throughfall  parameter  is  replaced  with  the  closely  related 
expression  1  -  Cf  in  order  to  reduce  the  number  of  parameters; 

2)  the  evaporation  factor  C/S  is  multiplied  by  Of  to  account 
for  the  asymptotic  limit  that  canopy  evaporation  vanishes  as 
the  canopy  vanishes;  and 

3)  n  is  chosen  to  be  less  than  unity  to  correspond  to  a  finite 
time  for  the  canopy  to  dry  following  rainfall  as  modelled 
in  Deardorff  (1978)  . 

Based  on  the  work  of  Leyton  et  al.  (1967),  a  value  of  n  *  1/2  is  inferred  which 
is  somewhat  less  than  the  n  =  2/3  value  chosen  by  Deardorff  (1978) . 

Once  the  canopy  is  saturated  (C»S),  all  additional  rainfall  is  assumed  to 
fall  through  to  the  ground.  This  is  analogous  to  assuming  that  drip  processes 
occur  instantaneously  so  that  the  canopy  is  never  temporarily  "supersaturated." 
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V.  A  Boundary  Layer  Formulation  for  Atmospheric  Models 


Abstract 


A  simple  model  for  the  atmospheric  boundary  layer  is  proposed  for  use  in 
operational  global  weather  prediction  models  and  other  models  where  simplicity 
is  required.  Surface  fluxes  are  represented  in  terms  of  similarity  theory  while 
turbulent  diffusivities  above  the  surface  layer  are  formulated  in  terms  of  bulk 
similarity  considerations  and  matching  conditions  at  the  top  of  the  surface 
layer.  The  boundary  layer  depth  is  represented  in  terms  of  a  modified  bulk 
Richardson  number.  Attention  is  devoted  to  the  interrelationship  between 
predicted  boundary  layer  growth,  the  turbulent  diffusivity  profile,  "counter- 
gradient"  heat  flux  and  truncation  errors.  The  model  is  especially  suited  for 
use  in  models  where  some  resolution  is  possible  within  the  boundary  layer,  but 
where  the  resolution  is  still  insufficient  for  resolving  the  detailed  boundary 
layer  structure. 

As  an  example  application,  the  model  is  used  to  study  the  influence  of 
surface  moisture  flux  on  the  boundary  layer  development. 
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t .  Introduction 

The  present  interest  in  medium  range  weather  prediction,  climate  response 
studies,  and  limited  area  forecast  models  creates  a  need  for  development  of 
single  formulations  of  surface  fluxes  and  transport  in  the  atmospheric  boundary 
layer.  A  number  of  boundary  layer  models  with  low  resolution  have  been  proposed 
for  use  in  large  scale  models.  One  approach  is  to  model  the  bulk  effect  of  the 
boundary  layer  by  interpolating  pertinent  variables  from  the  large  scale  model 
without  attempting  to  resolve  any  boundary  layer  structure  explicitly  (Clarke, 
1970;  Deardorff,  1972;  Smeda,  1979;  Chang,  1981). 

In  models  where  some  grid  levels  are  available  to  resolve  boundary  layer 
structure,  a  more  direct  approach  is  usually  adopted  by  expressing  turbulent 


diffusivities  in  terms  of  local  gradients  of  the  mean  profiles.  Models  of  this 
kind  have  been  used  mostly  in  cases  where  comparatively  high  resolution  is 
available  (column  models),  then  diffusivities  are  related  directly  to  the  local 
gradient  Richardson  number  (Zhang  and  Anthes,  1982),  or  to  a  Richardson  number 
adjustment  scheme  (Chang,  1979),  or  computed  in  conjunction  with  a  prescribed 
mixing  length  profile  (Busch  et  al . ,  1976;  Louis,  1979).  With  coarser  resolu¬ 
tion,  the  sensitivity  of  these  formulations  to  small  changes  in  the  mean  pro¬ 
files  becomes  a  disadvantage.  Inclusion  of  transport  terms  by  employing  the 
turbulence  energy  equation  (Mailhdt  and  Benoit,  1982;  Troen,  1982)  or  even 
higher  order  closure  schemes  (Yamada  and  Mellor,  1975;  Andrd  et  al.,  1978)  is 
presently  not  practical  for  use  in  large  scale  models  because  of  the  large 


computational  requirements. 

Here  we  develop  a  model,  where  turbulent  diffusivities  have  a  prescribed 
profile  shape  as  a  function  of  z/h  and  scale  parameters  derived  from  similarity 
arguments  where  z  is  the  height  above  ground  and  h  is  the  boundary  layer  top. 
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This  approach  partially  alleviates  resolution  requirements  yet  is  more  flexible 
them  the  purely  "bulk"  models. 

Similar  approaches  for  the  simulation  of  the  heated  boundary  layer  have 
been  proposed  by  Pielke  and  Mahrer  (1975),  and  Yu  (1977).  The  present  model, 
however,  differs  from  these  approaches  both  with  respect  to  the  profile  formula¬ 
tions  and  the  way  the  boundary  layer  height  is  determined. 

In  Section  6,  the  model  is  used  to  study  interactions  between  surface 
evaporation,  mixing  and  boundary  layer  growth  by  arbitrarily  specifying  the 
relationship  of  actual  surface  evaporation  to  the  potential  evaporation.  The 
intention  is  to  identify  certain  physical  feedbacks  which  have  not  been 
previously  considered. 

2 .  The  model 

a.  The  surface  boundary  layer 

The  surface  layer  parameter i ration  scheme  devised  by  Louis  (1979)  is  used 
to  relate  surface  fluxes  of  heat,  momentum  and  water  vapor  to  the  values  of 
temperature,  the  wind  components  and  specific  humidity,  all  at  the  lowest  model 
level.  The  layer  between  the  surface  and  the  lowest  model  level  is  thus 
considered  to  be  in  equilibrium  obeying  surface  layer  similarity.  The  basic 
advantage  of  this  formulation  is  computational  efficiency,  since  the  formulation 
avoids  an  iterative  process  which  is  otherwise  necessary  when  employing  the 
original  expressions  given  by  Businger  et  al.  (1971)  for  the  usual  range  of 
atmospheric  stability;  however,  the  correct  behavior  of  such  formulations  is 
uncertain  in  the  cases  of  extreme  stability  or  instability. 

b.  The  boundary  layer  above  the  surface  layer 

Above  the  surface  layer,  the  diffusion  equation  is  assumed  to  adequately 
describe  the  effect  of  turbulent  mixing  in  the  boundary  layer  except  for  the 
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modification  due  to  a  "counter-gradient"  term.  Thus,  model  closure  simplifies 
essentially  to  the  determination  of  the  diffusivity  profiles. 

As  in  Brost  and  Wyngaard  (1978),  the  momentum  diffusivity  is  modelled 
according  to  the  format 

K  -  u*  Jcz  <j>  - 1  ( 1  -  *)P  (1) 

m  *  ra  h 

where  u  is  the  surface  friction  velocity.  It  is  the  von  Karman  constant,  <f>  is 
*  m 

the  nondimensional  shear,  and  h  is  again  the  boundary  layer  height.  Eq.  (1)  is 
consistent  with  surface  layer  similarity  where 

K  =»  u.  kz  <b  *  (z/L)  for  z  «  h.  (2) 

m  •  m 

For  stable  conditions  we  use  ^  from  Businger  et  al.  (1971)  given  as 

<j>  «  1  +  4.7  z/L  (3) 

m 

where  I*  is  the  Obuknov  length. 

For  z  >  L  we  obtain  the  following  asymptotic  expression: 

K  =  0.09  L  u.  (1  -  £)P  z  >  L  .  (4) 

in  *  n 

That  is,  for  the  stable  case  L  becomes  the  relevant  length  scale  and  uA  the 
relevant  velocity  scale  for  the  entire  boundary  layer.  The  boundary  layer  depth 
h  here  only  enters  as  the  height  at  which  the  turbulence  varnishes  and  does  not 
influence  the  boundary  layer  velocity  scale. 

For  unstable  conditions 


m 


i 

m 


<|»  (z/L)  «  (1  -  7  z/L)_1/3,  z  «  h  .  (5) 

m 

The  exponent  of  -1/3  is  chosen  to  ensure  the  free  convection  limit  for  z  »  L. 
With  the  coefficient  chosen  to  be  7,  the  difference  between  Bq.  (5)  and  the 
original  expression  given  by  Businger  et  al.  (1971)  as  derived  from  the  Kansas 
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data  differs  by  less  than  6%  over  the  range  of  the  original  data  (-z/L  <  2). 

Here  we  consider  the  surface  layer  to  extend  upward  to  z  =  eh. 

Above  the  surface  layer  for  the  unstable  case,  we  assume  that  the  velocity 
scale  is  constant  and  given  by 

Ws  =  “  {U*3  +  7el™*3>1/3  <«> 

then  ( 1 )  becomes 

K  =  w  h  k  £  (1  -  ~)p,  z  >  eh  .  (7) 

m  s  n  n 

For  h  »  -L,  the  velocity  scale  approaches 

w  «  w*  (7ek)1/3  (8) 

s 

where  w  =  w*  a'  h) ^3  is  the  convective  velocity  scale.  For  e  *  .1,  and 
*  T  O 
O 

h  »  -L,  wg  approaches  .65  w#. 

The  expression  for  the  mixed  layer  velocity  scale  (6)  can  be  compared  with 
the  velocity  scale  developed  by  Hojstrup  (1982)  from  the  Kansas  and  Minnesota 
experiments.  Hojstrup's  expressions  for  the  velocity  variances  at  z/h  »  e  - 
. 1 ,  reduce  to 

w  -  (o  2  +  a  2  +  o  2)1/2  -  2.26  u  (1  +  2.75  (O.t  £)2/3)  1/2  .  (9) 

Apart  from  a  constant  factor  of  proportionality,  Bqs.  (6)  and  (9)  differ  by  less 
them  16%  in  the  range  of  h/-L  between  0  and  5000,  in  which  range  w  changes 

3 

by  a  factor  of  13. 

The  factor  (1-z/h)^  ensures  that  the  turbulent  mixing  approaches  zero 
towards  the  top  of  the  boundary  layer.  In  the  convective  case,  this  would 
correspond  to  a  capping  inversion  with  unimportant  turbulent  mixing  in  the 
overlaying  stratified  free  atmosphere.  The  eddy  diffusivity  is  zero  above  the 
boundary  layer  top. 


3.  Determination  of  boundary  layer  height 

For  well-defined  diurnal  variations,  a  rate  equation  for  growth  of  the 
daytime  mixed  layer  is  often  used  with  a  separate  model  for  the  depth  of  the 
nocturnal  boundary  layer.  However,  on  a  day-to-day  basis  at  an  arbitrary 
location,  the  top  of  the  daytime  boundary  layer  and  capping  inversion  are  often 
not  well-defined.  The  depth  of  the  boundary  layer  may  be  constantly  adjusting 
to  changing  synoptic  control  in  which  case  the  boundary  layer  growth  equations 
and  boundary  layer  collapse  provisions  become  too  specialized.  Even  in 
situations  with  definable  mixed  layer  growth,  the  vertical  resolution  in  most 
atmospheric  models  is  inadequate  to  define  a  capping  inversion. 

Alternatively,  Busch  et  al.  (1976)  choose  the  boundary  layer  top  to  be  the 
lowest  model  layer  where  the  gradient  Richardson  number  exceeds  a  critical 
value.  This  method,  however,  still  requires  good  resolution,  and  even  with  high 
resolution  may  lead  to  large  unphysical  oscillations  of  h  due  to  the  sensitivity 
of  the  gradient  Richardson  number  to  small  changes  in  the  mean  profiles. 

To  be  consistent  with  the  bulk  approach  adopted  in  Section  2,  we  determine 
the  boundary  layer  top  to  be  the  height  where  the  bulk  Richardson  number  over 
the  whole  boundary  layer  becomes  critical,  viz.. 


where  6v(h)  and  9a  are  the  virtual  potential  temperatures  at  the  boundary 


layer  top  and  near  the  surface,  respectively.  The  bulk  Richardson  number  is 


frequently  used  to  model  the  depth  of  the  stable  boundary  layer.  In  the  heated 


boundary  layer,  the  predicted  boundary  layer  top  occurs  just  above  the 


well-mixed  region.  The  thickness  of  the  implied  entrainment  region,  between  the 
well-mixed  region  and  predicted  boundary  layer  top,  depends  on  the  free  flow 


stratification  and  the  value  chosen  for  the  critical  bulk  Richardson  number. 
However,  the  downward  heat  flux  due  to  entrainment  implied  by  this  model  also 
depends  on  the  profile  of  diffusivity  in  this  region  and  the  way  such  a 
diffusivity  interacts  with  the  h-formulation  and  truncation  errors.  In  the 
case  of  vanishing  wind  speed,  and  thus  varnishing  shear  generation  of  turbulence 
by  the  mean  wind,  relationship  (10)  reduces  to  a  purely  thermodynaumic  argument 
as  applied  in  Holzworth  (1964). 

The  choice  of  the  lower  temperature  9g  in  (10)  plays  am  importamt  role. 

Since  the  largest  and  most  energetic  scales  of  turbulent  motion  in  the 

convective  boundary  layer  are  thermals,  it  seems  more  correct  to  define  9g  as  a 

measure  of  temperature  of  the  thermals  in  the  lowest  part  of  the  boundary 

layer.  This  cam  be  estimated  from  the  relevamt  velocity  scale  w  corresponding 

s 

to  a  thermal  turnover  time  of  h/ws.  The  virtual  temperature  excess  near  the 
surface  is  then 


(wir^r) 

V  o 


(11) 


where  (w' 9' )  is  the  surface  virtual  kinematic  heat  flux  amd  C  is  a  coefficient 
v  o 

of  proportionality.  Use  of  this  temperature  excess  to  estimate  the  boundary 
layer  top  from  (10)  could  be  viewed  as  a  parcel  approximation  which  neglects  the 
influence  of  entrainment  and  pressure  effects  on  the  thermal  ascent.  This 
overestimation  of  thermal  ascent  would  be  partially  compensated  by  neglect  of 
penetration  of  thermals  beyond  the  buoyancy  equilibrium  level.  The  attempt  to 
include  such  complexities  in  a  limited  resolution  model  would  not  be  appropriate 


due  to  large  truncation  errors 


For  simplicity,  the  temperature  excess  (11)  is  assumed  to  occur  at  the 


lowest  atmospheric  level  in  the  model,  z^ .  Then 


0=0  (z,)  +  0m  . 
S  v  1  T 


For  the  unstable  case,  the  modelled  boundary  layer  depth  (10-12)  depends  mainly 
on  the  thermodynamics  and  is  insensitive  to  the  choice  of  the  critical 
Richardson  number. 

However,  as  we  approach  neutral  conditions  from  the  unstable  side,  0T  in 
(11)  vanishes .  Then  0g  becomes  equal  to  the  virtual  temperature  at  the  first 
grid  level  and  the  modelled  boundary  layer  depth  depends  significantly  on  wind 


speed. 


Within  the  resolution  of  the  model,  the  boundary  layer  depth  relations 


(10-11)  compared  favorably  with  more  sophisticated  entrainment  models  for  the 


data  examined  here.  That  is,  for  sole  purpose  of  approximating  the  boundary 
layer  growth  rate,  the  main  role  of  the  convective  heating  can  be  captured 
rather  simply. 

Relationships  (11)  and  (12)  sure  not  relevant  for  the  stable  boundary 
layer.  Since  the  first  model  level  may  be  above  the  nocturnal  boundary  layer, 

0g  is  defined  to  be  the  surface  virtual  temperature  for  the  stable  case.  Here 
turbulence  is  generated  by  shear  only  and  the  Richardson  number  would  seem  to  be 
a  relevant  number.  Unfortunately,  the  turbulence  in  the  nocturnal  boundary 
layer  is  often  intermittent  and  layered  so  that  the  turbulence  at  a  given  level 
may  be  more  related  to  the  local  Richardson  number  them  the  bulk  Richardson 


number 


However,  in  low-resolution  large  scale  models,  the  structure  of  the 
nocturnal  boundary  layer  cannot  be  resolved  so  that  only  bulk  formulations  can 
be  implemented.  The  use  of  the  bulk  Richardson  number  (10)  to  predict  the  top 
of  the  nocturnal  boundary  layer  has  been  tested  in  a  number  of  studies  (some  of 
which  are  surveyed  in  Mahrt,  1981a),  and  provides  a  smooth  transition  to  the 
unstable  case  in  the  present  model. 

In  such  studies  the  critical  bulk  Richardson  number  is  typically  chosen 
between  .3  and  1.  Such  values  are  sometimes  tested  against  the  depth  of  the 
nocturnal  inversion  which  may  be  considerably  thicker  than  the  depth  of  the 
layer  of  continuous  turbulence  (Mahrt  et  al.,  1979;  Andr6  and  Mahrt,  1982).  On 
the  other  hand,  we  must  recognize  that  model  fluxes  imply  both  time  and 
horizontal  averaging.  Such  averaging  would  then  include  transport  induced  by 
mesoscale  and  terrain  induced  circulations  which  seem  to  be  important  in  the 
nocturnal  boundary  layer  even  over  very  weak  slopes  (Wyngaard  et  al.,  1978; 
Mahrt,  1981b).  Furthermore,  fluxes  may  occ  locally  in  space  and  time  even 
though  the  Richardson  number  evaluated  from  averaged  variables  is  large.  These 
factors  suggest  choosing  a  larger  critical  Richardson  number  for  computing  the 
boundary  layer  depth.  Here  we  choose  a  value  of  unity  for  the  critical  bulk 
Richardson  number. 

In  the  stable  case,  the  modelled  boundary  layer  depth  is  less  significant 
than  in  the  unstable  case  as  is  evident  from  comparison  of  (4)  and  (7).  That 
is,  the  implied  length  scale  of  the  mixing  asymptotically  becomes  independent  of 
h  and  proportional  to  L.  In  the  unstable  case,  the  value  of  h  significantly 
influences  the  length  scale  above  the  surface  layer;  however,  the  value  of  h 
becomes  insensitive  to  the  bulk  Richardson  number. 


(>*  .  .4  ii.  I'. 


4.  The  dif fusivities  for  heat  and  water  vapor 

We  make  the  usual  assumption  of  equating  dif fusivities  for  heat  and  water 
vapor.  The  turbulence  Prandtl  number  Pr  =  K^/K^  under  unstable  conditions  is 
found  to  be  strongly  dependent  on  stability  in  the  surface  layer  (Businger  et 
al . ,  1971).  In  the  mixed  layer  above  the  surface  layer,  the  Prandtl  number  is 
not  very  well-defined  because  local  gradients  may  vanish  and  fluxes  become  more 
related  to  bulk  gradients.  Thermals  and  boundary  layer  scale  eddies  transport 
properties  according  to  bulk  gradients  which  may  be  much  larger  than,  or  of 
opposite  sign  from,  gradients  in  the  boundary  layer  interior. 

The  simplest  way  to  include  this  nonlocality  is  to  incorporate  a  "counter- 
gradient"  term  as  discussed  by  Priestly  and  Swinbank  (1947)  and  Deardorff 
(1966).  Then  the  heat  flux  becomes 

wnr  *  (||  -  y)  •  (i3) 

In  the  present  formulation 

(w'F'T 

y“C-wTT^  *  (14) 

s 

This  prescription  is  consistent  with  the  formulation  of  the  temperature  excess 
for  thermals,  (11).  A  similar  correction  procedure  was  suggested  by  Deardorff 
(1973)  and  used  in  the  model  by  Mailhdt  and  Benoit  (1982),  except  that  the  free 
convection  velocity  scale  wA  was  chosen  to  be  the  velocity  scale.  The  use 
of  the  velocity  scale  wg  is  more  consistent  with  the  formulation  of  the  eddy 
exchange  coefficient  and  includes  the  reduction  of  thermal  buoyancy  by  mechani¬ 
cal  mixing.  Since  the  counter-gradient  factor  is  a  correction  to  the  flux 
predicted  by  the  eddy  exchange  coefficient,  the  value  of  y  and  therefore  C 
depend  on  the  formulation  of  the  eddy  diffusivity. 

This  interrelationship  is  clearest  at  the  level  where  the  potential  tem¬ 
perature  gradient  vanishes  as  it  reverses  with  height  from  weakly  unstable  to 


weakly  stable.  At  this  level,  say  z  *  z*,  (13)  becomes 


w' 9'  *  Kh(z*)y 


or  using  (14) 


hw 


Thus  C  cannot  be  specified  independently  of  the  formulation  of  as  has  been 
done  in  previous  studies.  Substituting  in  the  formulation  for  K^,  the  expres¬ 
sion  for  C  becomes 


-1 


C  =■  [pr  k  —  (l-z*/h)P] 


Note  that  C  becomes  independent  of  the  choice  of  the  velocity  scale  ws. 

Choosing,  for  example,  p  «*  2,  Pr  -  2  and  z*/h  -  1/2,  one  obtains  the  usual 
value  of  C  »  10.  We  will  adopt  this  value  for  subsequent  computations  although 
the  results  in  this  study  were  found  to  be  insensitive  to  the  numerical  value  of 
C. 

It  seems  necessary  to  adopt  a  counter-gradient  correction  term  for  trans¬ 
port  of  moisture  or  any  scaler,  since  thermals  also  transport  according  to  the 
bulk  moisture  gradient  and  therefore  create  flux,  even  where  the  local  mean 
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gradient  vanishes.  We  then  define  the  counter-gradient  factor  for  moisture  by 
assuming  that  it  can  be  formulated  with  the  same  coefficient  C  as  for  heat 
transport  in  which  case 


(15) 


where  (w'q* ’)  is  the  surface  moisture  flux.  An  analogous  correction  for 
momentum  is  not  adopted.  Because  of  the  pressure  effects,  thermals  cannot 
efficiently  transport  momentum  over  large  distances  and  the  gradient  of  momentum 
often  remains  significant  throughout  the  mixed  layer. 


a.  Prandtl  Number 

Busch  et  al.  (1976)  assume  the  Prandtl  number  to  obey  surface  layer 
relations  for  the  entire  boundary  layer  while  Mailhdt  and  Benoit  (1982)  assume 
that  the  Prandtl  number  is  independent  of  height  with  a  value  computed  at  4  m. 
In  the  present  development  we  match  heat  and  momentum  fluxes  at  the  top  of  the 
surface  layer  so  that 

u*9*  -  -*h  <4£  -  "«> 
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Combining  these  two  relationships,  using  the  definitions  of  nondimensional 
gradients  (fa,  fa),  substituting  for  y  from  (14)  and  solving  for  the  Prandtl 
number,  we  obtain 


(17) 


where  z  is  the  level  where  matching  (16)  is  applied,  here  taken  as  .lh.  Lacking 
other  evidence,  we  assume  that  the  Prandtl  number  is  constant  above  .lh. 
Expression  (17)  has  the  advantage  that  it  is  well  behaved  for  vanishing  Obukhov 
length  L. 

b.  Comparison  with  Wyngaard  and  Brost 

It  is  instructive  to  compare  the  diffusivity  profile  used  here  with  the 
diffusivity  profile  which  can  be  derived  from  the  large  eddy  simulations  by 
Wyngaard  and  Brost  (1983).  They  derive  expressions  for  the  gradients  of  a 
scalar  for  the  case  of  vanishing  entrainment  flux  (their  Eq.  (33))  and  for  the 
case  of  vanishing  surface  flux  (their  Bq.  (39)).  Profile  functions  were 
determined  by  comparing  with  numerical  simulations  for  the  case  -z/L  »  64.  When 
both  fluxes  are  present  we  can  obtain  the  effective  diffusivity  from  their 
individual  relationships.  Assuming  the  flux  in  the  boundary  layer  to  vary 
linearly  with  height,  the  diffusivity  relationship  derived  from  Wyngaard  and 
Brost  becomes 

(1  -  (1  -  R)  *) 

K  =  w .  h  - - - -  (18) 

*  r  a  -  .  o..  tj)-™ 

where  R  =  (wrcr)h  /(w*c' )  and  c  is  the  transported  quantity. 

As  discussed  by  Wyngaard  and  Brost  (1983),  this  form  is  not  well  behaved 
since  the  transport  by  large  eddies  or  thermals  is  not  directly  related  to  the 
local  gradients.  The  simplest  way  to  effectively  include  this  transport  process 
is  to  again  adopt  a  gradient  correction  factor  y  such  that 


where  the  prime  notation  on  the  coefficient  C*  indicates  that  the  free 
convection  velocity  scale  must  be  used  instead  of  ws  in  order  to  be 
consistent  with  (18).  Then  using  the  gradient  and  flux  from  Wyngaard  and  Brost 
(1983),  the  diffusivity  satisfying  (19)  and  the  gradient  correction  y  (14) 
becomes 

(1  -  (1  -  R)  £) 

K  ■  w#h - 2”- 3/2 - 2 --372 -  * 

R  (1  -  £)  3/2  +  0.4  (£)  3/2  +  C' 
h  h 

The  diffusivity  profile  modified  to  include  the  counter-gradient  correction  (20) 
is  shown  in  Pig.  1  for  different  values  of  R,  the  ratio  of  the  entrainment  flux 
to  the  surface  flux  .  C1  has  been  assigned  the  value  of  10  as  in  previous 
studies.  Eq.  (20)  and  Fig.  1  show  that  the  addition  of  the  counter-gradient 
correction  makes  the  diffusivity  well-behaved  and  only  moderately  sensitive  to 
the  value  of  the  ratio  R. 

c.  The  exponent  p 

In  the  case  of  the  diffusivity  for  heat,  the  ratio  R  is  often  approximated 
as  -.2  (Tennekes,  1973).  The  profile  for  the  heat  diffusivity  based  on  (17) 
for  values  of  the  exponent  p  between  2  and  3,  exhibits  good  agreement  with 


the  corresponding  profile  from  Wyngaard  and  Brost  (1983)  for  R=-.2  except  the 
heat  diffusivity  is  somewhat  larger  in  the  middle  of  the  boundary  layer  (Pig. 
2). 


Near  the  top  of  the  boundary  layer,  behavior  of  the  diffusivity  determines 
the  importance  of  the  entrainment  flux.  The  exact  relation  between  the  heat 
diffusivity  near  the  top  and  the  flux  ratio  R  is,  however,  complicated  by  the 
way  in  which  the  boundary  layer  height  is  determined,  and  is  in  practice  also 
found  to  be  sensitive  to  truncation  errors  in  the  numerical  model.  For  this 
reason  we  use  the  exponent  p  as  an  adjustable  parameter  chosen  to  give  reason¬ 
able  values  for  the  flux  ratio  R  in  a  particular  model  with  a  particular  resolu¬ 
tion.  It  is  obvious  that  the  exponent  p  could  be  dependent  on  stability.  For 
example,  Brost  and  Wyngaard  (1978)  choose  1.5  for  the  stable  case.  For  simpli¬ 
city,  and  lack  of  observational  evidence,  we  presently  specify  p  to  be 
independent  of  stability  with  a  constant  value  of  2. 

5 .  Numerical  technique 

The  time  integration  of  the  change  of  mean  profiles  due  to  boundary  layer 
turbulence  reduces  to  the  solution  of  the  diffusion  equation 


i*  =  i_  K  i i*  .  Y) 

3t  3z  K  '  3z  yt 


(21) 


where  X  «  [u,v,8,q].  At  each  time  step,  the  mean  profiles  are  used  first  to 
calculate  drag  coefficients  c^  and  c^  from  Louis  (1979)  and  then  calculate  the 
surface  fluxes.  Next,  boundary  layer  height  is  estimated  for  computation  using 
(10)  with  6S  *  0O  where  0O  is  the  surface  virtual  temperature  obtained 


from  surface  energy  balance  (Appendix).  Under  convective  conditions  (positive 
surface  heat  flux),  Eq.  (11)  is  then  used  together  with  (12)  and  (10)  to  obtain 
an  improved  estimate  for  the  boundary  layer  height  h.  Eq.  (21)  is  then  used  to 


step  forward  one  time  step  using  the  diffusivity  profile  (1).  The  layer  between 
the  ground  and  the  first  atmospheric  level  is  treated  as  an  equilibrium  surface 
layer  and  (21)  is  employed  only  above  this  layer.  The  boundary  conditions  thus 


become  flux  conditions  at  the  lowest  model  level.  The  boundary  layer  height,  h, 
counter-gradient  correction  term,  surface  values  of  temperature,  90,  and 
humidity,  <^j,  are  assumed  constant  during  each  time  step.  In  order  to  ensure 
computational  stability  with  large  time  steps,  it  is  necessary  to  use  an 
implicit  integration  technique;  here  we  use  the  fully  implicit  Cranck-Nicholson 
scheme  (or  Laasonen  scheme)  given  by 


(n+1)  (n)  _  A.  9  _  3  v(n+l) 
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where  the  superscripts  designate  the  time  level.  For  the  spatial  operator,  we 
employ  a  variation  of  the  finite-element  method  using  Chapeau  functions  for  the 
mean  variables .  Some  care  must  be  taken  to  ensure  that  no  unnecessary 
truncation  errors  are  introduced.  The  model  is  intended  to  be  used  with  only  a 
few  grid  levels  in  the  boundary  layer,  say  levels  at  50,  250,  500,  1000,  2000m. 
Therefore,  the  large  variation  of  the  diffusivities  near  the  top  of  the  boundary 
layer  will  typically  not  be  adequately  resolved.  To  minimize  truncation  errors 
from  finite  differencing  of  the  diffusivity  profile,  we  use  the  finite  element 
technique  with  the  analytical  form  of  K  from  (1).  The  method  can  be  developed 
as  follows:  suppose  X  is  written  in  terms  of  some  expansion  in  basic  functions 
ct^fz).  Multiplying  by  a^[z)  on  both  sides  of  (22)  and  integrating  from  the 
lowest  level  to  the  boundary  layer  top  yields 
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We  use  Chapeau  functions  for  the  basic  functions  defined  on  the  grid  as 


2  -  Vi 
zi  -  Vi 


for  *  <  z  <_  z^ 
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and  zero  outside  the  interval  -  z^+^ .  These  functions  span  all  piecewise 

linear  functions  on  the  grid.  Using  (24)  and  (1),  and  are  easily 
computed.  Solution  of  the  diffusion  step  then  becomes 

(\i  -  V1  xin+1  -  \i  xi”  •  >25> 
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The  expansion  coefficients  are  simply  the  gridpoint  values  by  virtue  of 


our  choice  for  in  (24). 


We  find  by  integration 
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All  other  terms  are  zero. 
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The  integrals  of  K  over  one  grid  interval  are  found  to  be  sufficiently  accurate 
when  approximated  by  the  center  value  multiplied  by  Az,  except  for  the  integral 
over  the  interval  containing  the  top  of  the  boundary  layer,  where  we  instead  use 
the  approximation 
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6.  Sensitivity  to  surface  moisture  flux 

One  of  the  goals  of  this  study  is  to  develop  a  boundary  layer  model  suffi¬ 
ciently  simple  to  allow  physical  analysis  of  the  qualitative  response  of  the 
boundary  layer  to  external  forcing.  In  this  section  we  will  attempt  to  identify 
interaction  between  boundary  layer  growth  and  surface  evaporation,  this  inter¬ 
action  is  normally  removed  in  mixed  layer  growth  studies  by  specifying  the 
surface  heat  flux. 

The  primary  influence  of  surface  evaporation  is  to  decrease  surface  heating 
which  in  turn  reduces  the  boundary  layer  growth  rate  from  that  of  the  case  of  no 
evaporation.  However,  this  interaction  is  limited  by  negative  feedback  since 
the  surface  evaporation  also  acts  to  increase  the  boundary  layer  moisture. 
Furthermore,  reduced  surface  heating  resulting  from  surface  evaporation  leads  to 
less  turbulence,  less  downward  mixing  of  dryer  air  and  lower  surface  wind  speed 
due  to  weaker  downward  mixing  of  momentum.  These  influences  all  lead  to  smaller 
potential  evaporation  compared  to  the  case  of  a  dry  surface. 

In  order  to  speculate  on  the  relative  importance  of  these  influences,  we 
now  examine  the  response  of  boundary  layer  development  to  the  surface  evapora¬ 
tion  using  the  model  constructed  above.  We  specify  the  initial  conditions  to  be 
the  observed  atmosphere  at  0600  Local  Standard  Time  on  Wangara  day  33  (Clarke  et 
al . ,  1971).  These  are  used  as  arbitrary  initial  conditions  since  the  lack  of 
moisture  flux  data  preclude  actual  comparisons  with  the  Wangara  data.  The  model 
is  integrated  for  a  48-hour  period  using  the  net  surface  radiation  reported  in 
Clarke  et  al.  (1971).  The  vertical  resolution  is  50m  with  a  vertical  domain  of 
2  km  adequate  for  present  purposes.  Varying  the  resolution  by  a  factor  of  two 
exerted  little  influence  on  the  results. 


As  the  specified  water  availability  factor  $  =  E/Bp  increases  from  zero 
(no  surface  evaporation)  to  .25  and  then  .5,  the  daytime  boundary  layer  growth 
rate  decreases  significantly  (Figs.  3-4).  Further  increases  of  (3  lead  to 
progressively  smaller  changes.  This  nonlinear  dependence  of  boundary  layer 
growth  rate  on  surface  moisture  availability  is  due  to  the  reduction  of 
potential  evaporation  (Fig.  4)  discussed  above.  Of  importance  is  that  the 
difference  between  the  present  simple  model  of  boundary  layer  depth  and  more 
sophisticated  mixed  layer  growth  models  is  small  compared  to  the  influence  of 
possible  variations  of  the  influence  of  surface  evaporation. 

Decreasing  the  atmospheric  moisture  increases  the  potential  evaporation 
which  in  turn  reduces  the  surface  heating  and  boundary  layer  growth  rate. 
Decreasing  the  specific  humidity  by  a  factor  of  two  exerts  only  a  modest  influ¬ 
ence  on  the  boundary  layer  growth  rate  for  the  Wangara  conditions  (Fig.  5). 

Using  the  boundary  layer  growth  rate  as  an  indicator  of  boundary  layer 
development.  Figs.  6-7  show  the  expected  influence  of  variations  of  atmospheric 
stratification  and  solar  radiation.  Varying  the  coefficient  C  for  the  gradient 
correction  (11)  by  a  factor  of  two  caused  negligible  modification  of  the  growth 
rate. 

7.  Cloud  scheme 

The  boundary  layer  model  includes  a  cloud  scheme  based  on  the  model  at  the 

European  Center  for  Medium-range  Heather  Forecasting  (Geleyn,  1980).  The 

fractional  boundary  layer  cloud  cover  is  determined  as  a  function  of  the  excess 

of  the  relative  humidity  over  a  critical  value  (here  .5)  such  that 
CC  -  [2  (RH  -  .5)  ]  RH  _>  .5 

0  RH  <  .5 

where  CC  is  the  fractional  cloud  cover  and  RH  is  the  relative  humidity. 


h 


TIME  (HRS) 


Fig .  3 . 


The  boundary  layer  depth  as  a  function  of  local  standard  time  for 
E/Ep=0  (open  circles),  .25  (open  squares),  .5  (solid  triangles), 
.75  (open  triangles)  and  1  (solid  circles). 
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Fig.  4.  Boundary  layer  depth  and  potential  evaporation  as  a  function  of 
specified  E/Ep  for  1200  local  time  for  day  33  (solid  line)  and 
day  34  (dashed  line ) . 


*lv 


:  _  a_: 


fixate 


TIME  (HRS) 


Fig.  5.  Boundary  layer  depth  as  a  function  of  local  time  for  the  Wangara 

initial  conditions  (open  circles)  and  for  the  case  with  the  initial 
specific  humidity  decreased  by  a  factor  of  two  (solid  circles), 
B/Ep-1. 


TIME  (HRS) 


Boundary  layer  depth  as  a  function  of  local  time  for  the  prototype 
experiment  (solid  circle),  solar  radiation  increased  by  50%  (open 
squares)  and  solar  radiation  reduced  by  a  factor  of  two  (open 
circles ) . 
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In  the  present  model  this  expression  is  evaluated  at  the  second  boundary  layer 
level  above  the  surface. 

He  have  not  objectively  tested  this  scheme  since  the  general  format  has 
already  been  subject  to  considerable  scrutiny.  Our  opinion  is  that  this  cloud 
representation  can  probably  be  improved  considerably,  although  needed  simplicity 
and  truncation  errors  will  limit  the  degree  of  potential  improvement.  The  most 
obvious  generalization  would  be  to  express  the  critical  relative  humidity  as  a 
function  of  basic  cloud  type. 

8 .  Conclusions 

This  study  develops  a  numerical  model  of  the  atmospheric  boundary  layer 
which  requires  only  modest  vertical  resolution  and  is  sufficiently  simple  for 
use  in  concert  with  other  models.  For  example,  the  formulation  of  the  boundary 
layer  depth  includes  the  main  features  of  mixed  layer  growth  but  does  not  lead 
to  special  interpolation  problems  occurring  with  use  of  sophisticated  mixed 
layer  growth  models  with  low  vertical  resolution.  The  simple  boundary-layer 
depth  model  adopted  here  includes  the  influences  of  mixing  generated  by  both 
shear  and  surface  heating  and  allows  for  a  smooth  transition  to  the  stable  case. 

The  "counter-gradient  correction"  to  the  heat  transport  by  boundary  layer 
scale  convective  eddies  has  been  generalized  to  be  consistent  with  surface  layer 
similarity  theory  and  at  the  same  time  allows  for  continuous  transition  to  the 
mechanically  mixed  boundary  layer. 

As  an  example  application,  the  model  is  used  to  study  the  response  of  the 
boundary  layer  to  the  surface  moisture  flux.  This  flux  is  formulated  as  a 
specified  fraction  of  the  potential  evaporation.  The  potential  evaporation  is 
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represented  by  a  Penman  expression.  The  daytime  growth  of  the  modelled  boundary 
layer  is  reduced  significantly  by  modest  availability  of  surface  moisture. 
However,  additional  surface  moisture  leads  to  less  important  changes  since  the 
return  influence  of  the  boundary  layer  modifications  leads  to  reduced  potential 
evaporation.  These  tentative  results  of  the  surface-atmosphere  coupling  suggest 
that  distinction  between  modest  surface  moisture  deficit  and  large  surface 
moisture  deficit  is  more  important  than  distinction  between  a  modest  moisture 
deficit  and  saturated  surface  conditions.  This  sensitivity  to  substantial 
moisture  deficits  would  be  augmented  by  inclusion  of  the  nonlinear  dependence  of 
transpiration  on  the  soil  moisture  deficit. 


Appendix 


We  assume  that  the  model  input  from  the  radiation  calculations  in  the 
larger  scale  model  includes  both  short  and  longwave  contributions.  The  surface 
energy  balance  then  becomes 

(1  -  a)  S+  +  L+  -  OT4  =  G  +  H  +  LE  (Al) 

o  o 

where  S+  is  the  downward  shortwave  radiation,  a  is  the  albedo,  L+  downward 
longwave  radiation,  G  the  ground  heat  flux  positive  when  downward,  Hq  and  LE 
are  the  sensible  and  latent  heat  fluxes  respectively,  positive  when  upward. 

The  latent  heating  is  related  to  the  potential  evaporation  through  the 
relationship 

EH0  Bp 

where  is  the  potential  evaporation  rate.  The  coefficient  3  is  related  to 
the  soil  moisture  deficit  and  plant  resistance  to  transpiration  (Monteith, 
1981).  The  present  discussion  does  not  explicitly  represent  the  contributions 
to  $.  Instead  8  is  arbitrarily  varied  in  order  to  study  the  response  of  the 
atmospheric  boundary  layer  to  the  surface  moisture  flux. 

The  computation  of  Ep  follows  the  modified  Penman  method  (of  Chapter  II) 
except  for  one  modification  motivated  by  numerical  efficiency.  The  original 
Penman  method  requires  knowledge  of  the  net  radiation  and  eliminates  surface 
temperature  as  a  parameter.  Here,  however,  we  need  the  surface  temperature  for 
radiation  and  surface  heat  flux  calculations,  and  we  therefore  use  (Al)  to 
determine  T0  in  the  manner  outlined  below  (A2-A3). 

In  Chapter  II,  it  was  necessary  to  avoid  the  use  of  surface  temperature 
since  that  development  defines  potential  evaporation  for  comparison  with  obser- 
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vations  of  evapotranspiration  carried  out  in  Chapter  IV.  In  such  cases,  use  of 
surface  temperature  over  land  must  be  avoided  since  it  cannot  be  defined  in  a 


useful  observational  framework.  Thus  the  usual  combination  of  the  surface 
energy  balance  with  the  bulk  aerodynamic  relationship  was  performed  in  Chapter 
II  in  order  to  eliminate  surface  temperature. 

Once  the  empirical  formulations  were  completed  in  Chapter  IV,  the  surface 
energy  balance  and  surface  temperature  are  reintroduced  in  the  modelling 
situation  of  this  chapter.  It  is  not  inconsistent  to  use  the  Penman  formulation 


in  concert  with  the  surface  energy  balance,  since  the  mathematical  situation  of 
two  unknowns  and  two  equations,  is  preserved.  The  usual  Penman  method  is  based 


on  a  linearization  of  the  saturation  vapor  pressure  in  terms  of  the  temperature 
difference  between  the  lowest  atmospheric  level  and  the  surface.  Linearizing 
the  upward  radiation  we  obtain 

°ro4  *  (1  +  <  (  ° ■T-1))  (A2) 


and  using  the  Penman  method  leads  to  the  following  expression  for  the  potential 
evaporation: 


K  h 

E  =  [a(1  +  r)  +  RA]  (A  +  1  +  r)' 

P  * 
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=  exchange  coefficient  for  heat 


p  =  density 


L  =  heat  of  evaporation  for  water 


•  >  •> 


0^  =>  specific  heat  of  air  at  constant  pressure 
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0  =  Stefan-Bolzman  constant 


p  =  surface  pressure 


Rgas  =  9as  constant  assumed  equal  to  dry-air  value 
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e  «  ratio  of  water  molecular  weight  to  that  of  dry  air  =  0.622 


q,  *  specific  humidity  at  atmospheric  level 


q^  -  saturation  specific  humidity  at  temperature 


T,  ■  temperature  at  atmospheric  level 
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(Al)  is  then  used  to  estimate  T 


T  =  T,  + 
o  1 


[r(.]_I_B  +  1)  -  &a](A  +  1  +  r) 


-1 


(A4) 


As  in  previous  models,  we  do  not  distinguish  between  the  surface  air 
temperature  at  the  level  of  the  roughness  elements  as  used  in  surface  layer 
similarity  theory  and  the  effective  surface  radiation  temperature.  Without  this 
assumption,  we  would  need  a  more  detailed  formulation  of  surface  conditions  than 
is  included  in  the  present  development.  It  should  also  be  noted  that  model 
surface  temperatures  have  an  uncertain  relationship  to  the  actual  temperature  of 
the  soil  surface  or  plant  canopy. 
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VI .  Concluding  Remarks 

The  research  described  in  this  report  develops  a  model  for  stability- 
dependent  potential  evaporation  (Chapter  II),  soil  hydrology  (Chapter  III), 
transpiration  and  the  canopy  water  budget  (Chapter  IV)  and  the  atmospheric 
boundary  layer  (Chapter  V) .  The  construction  of  these  models  was  partly 
motivated  by  the  important  influences  of  the  surface  moisture  flux  on  the 
diurnal  development  of  the  boundary  layer,  related  cloud  fields  and  flow 
patterns  in  general.  Since  existing  global  circulation  models  neglect  the 
influence  of  surface  moisture  flux  over  land  or  incorporate  them  in  a  cursory 
manner,  much  of  the  model  development  in  this  report  is  original  in  nature. 

The  software  logistics  of  the  combination  of  these  models  is  described  in  the 
User's  Manual  corresponding  to  this  report.  A  rough  flow  chart  is  provided  in 
Fig.  1. 

The  physics  of  the  interaction  between  models  is  primarily  tied  to  the 
hydrological  budget.  Surface  moisture  flux  from  the  soil  and  plant  canopy  leads 
to  a  more  moist  and  cooler  boundary  layer.  The  difference  between  a  moist  and 
dry  surface  can  cause  major  differences  in  the  diurnal  evolution  of  the  atmos¬ 
pheric  boundary  layer  over  land  as  was  documented  with  the  boundary  layer  model 
developed  in  Chapter  V. 

The  modification  of  the  boundary  layer,  in  turn,  influences  the  surface 
moisture  flux.  The  atmospheric  humidity  deficit,  wind  speed  and  transmission  of 
shortwave  radiation  control  the  "maximum  allowed"  surface  moisture  flux.  This 
control  is  usually  formulated  through  the  Penman  expression  for  potential 
evaporation.  Chapter  II  develops  a  stability-dependent  expression  for  potential 
evaporation  for  use  in  atmospheric  models.  Previous  models  of  potential 
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solar  radiation 
wind 

humidity  < - turbulent 


FLOW  CHART  FOR  MOISTURE 


Fig.  la.  Diagram  of  terms  in  moisture  budget  (double  arrows)  and 
interaction  with  turbulent  boundary  layer  ( thin  arrows ) . 


evaporation  are  either  too  complicated  for  use  in  general  circulation  models  or 
have  omitted  the  important  influence  of  stability  on  potential  evaporation. 


! 

As  noted  in  Chapter  I ,  the  total  surface  evaporation  can  be  written  as 
E  =»  +  E>j>  +  Eq  ■  F  Ep 

where  %ir  is  again  the  direct  evaporation  from  the  soil,  Erp  is  the  trans¬ 
piration,  Ec  is  the  evaporation  of  precipitation  water  intercepted  by  the 
canopy,  Ep  is  the  potential  evaporation  and  P  is  a  function  which  relates  the 
various  contributions  of  the  surface  moisture  flux  to  the  potential 
evaporation . 

The  function  P  can  be  written  as 


z  — z 
2  1 


+  Vf^V^  g(0l)  +  g( ®2 )  1  [1-<C/S)n]  +  Of(C/S)n 

0  -  0  ... 

=  w.lt 

9(0)  *  0  -  0  ... 

ref  wilt 


The  function  Psoil  is  not  formally  defined  but  rather  represents  the  modelled 
process  whereby  the  direct  evaporation  from  the  soil  is  at  the  potential  rate 
until  the  corresponding  surface  soil  moisture  decreases  below  the  air-dry  value 
as  discussed  in  Chapter  III,  Section  4.  At  this  point,  the  soil  evaporation 
rate  becomes  less  than  the  potential  rate  and  becomes  controlled  by  the  soil 
moisture  profile  and  surface  air-dry  value.  The  relevant  soil  moisture  profile 
is  determined  by  the  soil  moisture  in  a  thin  upper  layer  of  the  soil  of  depth 
.  Transpiration  is  more  dominated  by  the  soil  moisture  in  a  deeper  layer  of 
depth  Z2~*l« 


The  contributions  to  the  surface  moisture  flux  for  the  total  grid  area  are 
modulated  by  the  fractional  coverage  of  vegetation.  Of,  which  reduces  direct 
evaporation  from  the  soil  and  leads  to  non-zero  transpiration  or  evaporation  of 
precipitation  intercepted  by  the  canopy.  The  intercepted  water,  C,  in  turn, 
reduces  the  transpiration  according  to  the  amount  of  canopy  water  with  respect 
to  the  total  water  capacity  of  the  canopy,  S. 

The  transpiration  is  controlled  by  the  stomatal  resistance  of  the  vegeta¬ 
tion  represented  by  the  plant  coefficient,  ky,  and  controlled  by  the  soil 
moisture  deficit  represented  by  the  function  g(0).  Transpiration  becomes  mono- 
tonically  reduced  after  the  soil  volumetric  water  content,  9,  decreases  below 
the  reference  value,  9ref.  Transpiration  vanishes  altogether  when  the  soil 
water  content  is  depleted  to  the  wilting  value,  9wilt*  Again,  transpiration 
is  controlled  mainly  by  the  water  content  of  the  deep  lower  layer,  02 ,  while 
direct  evaporation  from  the  soil  is  dictated  by  the  soil  water  content  of  the 
thin  upper  layer,  0j_. 

Numerical  experiments  were  performed  where  the  soil,  canopy  and  atmospheric 
boundary  layer  models  were  run  in  concert.  These  sensitivity  tests  are  not 
reported  here.  Significant  meaning  could  not  be  attached  to  the  results  except 
that  the  total  model  appeared  to  behave  in  a  reasonable  manner.  Attempts  to 
compare  the  model  with  classical  soil  moisture  data  from  the  USDA  laboratory  in 
Phoenix  led  to  identification  of  severe  inconsistencies  in  the  corresponding 
atmospheric  data  needed  to  evaluate  the  model.  Comparisons  with  data  from  the 
o'Niell  field  program  were  also  conducted;  however,  data  gaps  and  the  limited 
availability  of  the  soil  moisture  data  became  sufficiently  important  to  render 
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comparisons  inconclusive.  Recent  data  from  the  KMNI  facility  near  de  Bilt, 
Netherlands  may  provide  more  meaningful  comparisons  although  this  data  repre¬ 
sents  a  rather  narrow  range  of  soil  moisture  conditions. 

Certainly  considerable  effort  is  required  to  further  examine  the  inter¬ 
action  between  the  various  model  components  for  a  variety  of  atmospheric  condi 
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A  method  of  calculating  surface  evapotranspiration  by  separately  including 
the  effects  of  vegetation  and  atmospheric  evaporative  demand  under  the  condition 
of  nonlimiting  soil  moisture  is  presented.  A  literature  survey  is  conducted  to 
determine  the  effects  of  plants  on  evapotranspiration. 

To  represent  the  atmospheric  evaporative  demand,  the  original  potential 
evaporation  equation  of  Penman  (1948)  is  utilized  and  then  modified  to  include 
the  effect  of  atmospheric  stability  using  turbulent  exchange  coefficients 
formulated  by  Louis  et  al .  (1982).  The  original  and  modified  Penman  expressions 
are  compared  for  different  asymptotic  cases.  Using  boundary  layer  data  from  the 
Wangara  experiment  (Clarke  et  al.,  1971),  the  diurnal  variations  of  the  original 
and  modified  Penman  equations  are  compared.  The  daily  total  potential 
evaporation  using  linearized  and  integrated  forms  of  the  original  and  modified 
expressions  are  also  compared.  Finally,  the  nonlinear  effects  of  averaging  both 
the  original  and  modified  expressions  sure  examined.  It  is  found  that  including 
the  diurnal  variations  of  stability  in  the  modified  expression  causes  large 
hourly  differences  with  the  original  expression  under  non-neutral  conditions, 
while  daily  averages  of  the  two  compared  fairly  well.  The  diurnal  variation  of 
the  surface  moisture  flux  appears  to  be  much  larger  than  predicted  by  the 
original  Penman  expression.  However,  the  original  Penman  expression  remains  a 
reasonable  estimate  of  the  24-hour  total  potential  evaporation. 
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Abstract 


The  Penman  relationship  for  potential  evaporation  is  modified  to  simply 
include  the  influence  of  atmospheric  stability  on  turbulent  transport  of  water 
vapor.  Explicit  expressions  for  the  stability-dependent,  surface  exchange 
coefficient  developed  by  Louis  are  used.  The  diurnal  variation  of  potential 
evaporation  is  computed  for  the  stability-dependent  and  original  Penman 
relationship  using  Wangara  data. 

The  influence  of  afternoon  instability  increases  the  aerodynamic  term  of  the 
modified  Penman  relationship  by  50%  or  more  on  days  with  moderate  instability. 
However,  the  unmodified  Penman  relationship  predicts  values  of  daily  potential 
evaporation  close  to  that  of  the  stability-dependent  relationship.  This 
agreement  is  partly  due  to  compensating  overestimation  during  nighttime  hours. 
Errors  due  to  use  of  daily-averaged  variables  are  examined  in  detail  by 
evaluating  the  nonlinear  interactions  between  the  diurnal  variation  of  the 
variables  in  the  Penman  relationship. 

A  simpler  method  for  estimating  the  exchange  coefficient  is  constructed  from 
an  empirical  relationship  between  the  radiation  Richardson  number  and  the 
Obukhov  length.  This  method  is  less  accurate,  but  allows  estimation  of  the 
stability-dependent  exchange  coefficient  using  only  parameters  already  required 
for  evaluation  of  the  Penman  relationship.  Finally,  the  diurnal  variation  of 
the  atmospheric  resistance  coefficient  appearing  in  the  Penman-Monteith 
relationship  is  presented. 
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Abstract 


A  two-layer  model  of  soil  hydrology  is  developed  for  applications  where  only 
limited  computer  time  and  complexity  are  allowed.  Volumetric  soil  water  is 
computed  in  a  thin  upper  layer  for  use  in  calculation  of  surface  evaporation. 
Storage  of  water  is  computed  for  an  underlying  deeper  layer. 

In  an  effort  to  identify  the  influence  of  significant  asymmetric  truncation 
errors  in  the  two-layer  model,  this  model  is  compared  with  the  100-level  model 
of  Boersma  et  al.  (1983).  Comparisons  are  made  for  modelled  soils  with  clay, 
loam  and  sand  properties  for  various  time  dependencies  of  potential  evaporation 
and  precipitation.  Truncation  errors  in  the  resulting  two-layer  model  appear  to 
be  modest  at  least  compared  to  errors  associated  with  difficulty  in  estimation 
of  the  hydraulic  diffusivity  and  its  strong  dependence  on  soil  water  content. 

Minimization  of  the  influence  of  truncation  errors  requires:  1)  choosing 
the  upper  layer  to  be  sufficiently  thin,  2)  allowing  the  soil  water  gradient  to 
directly  control  surface  evaporation  and  3)  using  suitable  numerical 
implementation  of  the  evaluation  of  internal  soil  water  flux. 
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ABSTRACT 


A  simple  model  for  the  atmospheric  boundary  layer  is  proposed  for  use  in 
operational  global  weather  prediction  models  and  other  models  where  simplicity 
is  required.  Surface  fluxes  are  represented  in  terms  of  similarity  theory  while 
turbulent  dif fusivities  above  the  surface  layer  are  formulated  in  terms  of  bulk 
similarity  considerations  and  matching  conditions  at  the  top  of  the  surface 
layer.  The  boundary  layer  depth  is  represented  in  terms  of  a  modified  bulk 
Richardson  number.  Attention  is  devoted  to  the  interrelationship  between 
predicted  boundary  layer  growth,  the  turbulent  diffusivity  profile,  "counter¬ 
gradient"  heat  flux  and  truncation  errors.  The  model  is  especially  suited  for 
use  in  models  where  some  resolution  is  possible  within  the  boundary  layer,  but 
where  the  resolution  is  still  insufficient  for  resolving  the  detailed  boundary 
layer  structure. 

As  an  example  application,  the  model  is  used  to  study  the  influence  of 


surface  moisture  flux  on  the  boundary  layer  development 
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